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Abstract 

Recently,  a  framework  for  multiscale  stochastic  modeling  was  introduced  based  on 
coarse-to-fme  scale-recursive  dynamics  defined  on  trees.  This  model  class  has  some 
attractive  characteristics  which  lead  to  extremely  efficient,  statistically  optimal  signal 
and  image  processing  algorithms.  In  this  paper,  we  show  that  this  model  class  is  also 
cjuite  rich.  In  particular,  we  describe  how  1-D  Markov  processes  and  2-D  Markov 
random  fields  (MRF’s)  can  be  represented  within  this  framework.  Markov  processes 
in  one-dimension  and  Markov  random  fields  in  two-dimensions  are  widely  used  classes 
of  models  for  analysis,  design  and  statistical  inference.  The  recursive  structure  of  1-D 
Markov  processes  makes  them  simple  to  analyze,  and  generally  leads  to  computationally 
efficient  algorithms  for  statistical  inference.  On  the  other  hand,  2-D  MRF’s  are  well 
known  to  be  very  difficult  to  analyze  due  to  their  non-causal  structure,  and  thus  their 
use  typically  leads  to  computationally  intensive  algorithms  for  smoothing  and  parameter 
identification.  Our  multiscale  representations  are  based  on  scale-rectirsive  models,  thus 
providing  a  framework  for  the  development  of  new  and  efficient  algorithms  for  Markov 
processes  and  MRF’s.  In  1-D,  the  representation  generalizes  the  mid-point  deflection 
construction  of  Brownian  motion.  In  2-D,  we  use  a  further  generalization  to  a  “mid-line” 
deflection  construction.  Our  exact  representations  of  2-D  MRF’s  are  of  potentially  high 
dimension,  and  this  motivates  a  class  of  approximate  models  based  on  one-dimensional 
wavelet  tran.sforms.  We  demonstrate  the  use  of  these  models  in  the  context  of  texture 
representation  and  in  particular,  we  show  how  they  can  be  used  as  approximations  for 
or  alternatives  to  well-known  MRF  texture  models.'  We  illustrate  how  the  quality  of 
the  representations  varies  as  a  function  of  the  underlying  MRF  and  the  complexity  of 
the  wavelet-based  approximate  representation. 
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DL-H-418524,  by  the  Office  of  Naval  Research  under  Grant  N00014-91-J-1004,  by  the  Army  Research  Office 
under  Grant  DAAL03-92-G-0115,  and  by  the  Air  Force  Office  of  Scientific  Research  under  Grant  AFOSR- 
92-J-0002.  The  last  two  authors  were  also  supported  by  the  AFOSR  under  contract  F49620-91-C-0047. 

^M.I.T  Laboratory  for  Information  and  Decision  Systems,  Cambridge,  MA  02139.  Electronic  mail  ad¬ 
dress;  luettgen@athena.mit.edii 
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1  Introduction 


In  this  paper,  we  describe  how  to  use  a  class  of  multiscale  stochastic  models  to  represent  1-D 
Markov  and  reciprocal  processes  and  2-D  Markov  random  fields  (MRF’s).  Markov  models  in 
one  dimension  provide  a  rich  framework  for  modeling  a  wide  variety  of  biological,  chemical, 
electrical,  mechanical  and  economic  phenomena  [10].  Moreover,  the  Markov  structure  makes 
the  models  very  simple  to  analyze,  so  that  they  often  can  be  easily  ai>plied  to  statistical 
inference  problems  (such  as  detection,  parameter  identification  and  state  estimation)  as  well 
as  problems  in  system  design  (e.g.  control  and  queuing  systems). 

In  two  dimensions,  MR.F’s  also  have  been  widely  used  as  models  for  physical  systems 
[.3,  4,  .51,  31],  and  more  recently  for  images.  For  example,  Gaussian  fields  [64]  have  been 
used  as  image  texture  models  [22,  23,  36,  13,  47,  48],  and  the  more  general  Gibb’s  fields 
have  been  used  as  prior  models  in  image  segmentation,  edge  detection  and  smoothing  prob¬ 
lems  [5,  33,  52,  49].  Causal  sub-classes  of  MRF’s,  such  as  Markov  Mesh  Random  Fields 
[1,  28]  and  Non-Symmetric  Half-Plane  Markov  chains  [35]  lead  to  two-dimensional  versions 
of  Kalman  filtering  algorithms  when  the  fields  axe  Gaussian  [65].  In  addition,  efficient  fast 
Fourier  transform  algorithms  are  available  for  stationary  Gaussian  fields  defined  on  toroidal 
lattices  [22,  36,  12].  In  general,  however,  Markov  random  field  models  lead  to  computa¬ 
tionally  intensive  algorithxns  (e.g.  stochastic  relaxation  [33])  for  estimation  problems.  In 
addition,  parameter  identification  is  difficult  for  MRF  models  due  to  the  problem  of  com¬ 
puting  the  partition  function  [4,  37,  53].  Thus,  while  Markov  random  fields  provide  a  rich 
structure  for  multidimensional  modeling,  they  do  not  generally  lead  to  the  simple  analysis 
and  computationally  efficient  algorithms  that  1-D  Markov  pi-ocesses  do. 

These  computational  issues  are  the  most  important  obstacle  to  the  application  of  MRF 
models  to  a  broader  range  of  problems,  and  are  the  principal  motivations  for  the  in¬ 
vestigation  in  this  paper  of  the  richness  of  the  class  of  multiscale  stochastic  processes 
[16,  17,  18,  19,  8,  9],  and  in  particular  of  how  such  multiscale  processes  can  be  used  to 
exactly  and  approximately  represent  Markov  random  fields.  Our  multiscale  stochastic  pro¬ 
cesses  are  described  by  scale-recursive  models,  which  lead  naturally  to  computationally 
efficient  scale-recursive  algorithms  for  a  variety  of  estimation  problems.  For  instance,  fast 
smoothing  algorithms  are  developed  for  a  class  of  Gaussian  px-ocesses  by  Chou  et.  al.  in 
[16,  17,  18,  19].  Also,  Bouman  and  Shapiro  demonstrate  how  a  related  multiscale  discrete 
random  field  [8,  9]  leads  to  an  efficient  sequexitial  MAP  estimator.  In  this  paper,  we  dexnon- 
strate  how  a  simple  generalization  of  the  models  in  [16,  17,  18,  19]  leads  to  classes  of  models 
which  can  be  used  to  represent  all  1-D  Markov  px-ocesses  and  2-D  Markov  randoxn  fields. 
The  significance  of  this  x’esult  is  not  only  that  it  opexis  the  door  to  the  possibility  of  new  and 
efficient  algorithxns  for  MRF  models,  but  also  that  it  suggests  that  this  multiscale  modeling 
framework  may  be  a  decidedly  superior  basis  for  image  and  random  field  xuodeling  and 
analysis  than  MRF’s  both  because  of  the  efficient  algorithxns  it  admits  and  because  of  the 
rich  class  of  phenomena  it  can  be  used  to  describe. 

Our  multiscale  representations  of  1-D  recipx'ocal  processes  and  2-D  MRF’s  rely  on  a 
generalization  of  the  mid-point  deflection  technique  for  constructing  a  Brownian  motion  in 
one  dixnension  [27,  32,  41].  To  construct  a  Brownian  motion  sample  path  over  an  interval 
by  mid-point  deflection,  we  start  by  randomly  choosing  values  for  the  process  at  the  two 
boundary  points,  i.e.  the  initial  and  final  points,  of  the  interval  according  to  the  joint 
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Figure  1:  The  first  three  levels  of  a  “mid-point  deflection”  construction  of  a  Brownian  motion 
sample  path  are  shown  on  the  left.  The  construction  generates  a  sequence  of  approximations 
based  on  linear  interpolations  of  samples  of  the  Brownian  motion  at  the  dyadic  points.  On 
the  right,  the  basis  functions,  integrals  of  the  Haar  wavelet,  in  this  construction  are  shown. 

probability  distribution  implied  by  the  Brownian  motion  model.  We  then  use  these  two 
values  to  compute  the  expected  value  of  the  process  at  the  mid-point,  which  is  just  the 
average,  and  then  add  to  that  a  Gaussian  random  variable  with  zero  mean  and  variance 
equal  to  the  variance  of  the  error  in  this  mid-point  prediction.  These  steps  are  illustrated  in 
the  left  side  of  Figure  1.  The  process  is  then  continued  by  using  the  original  boundary  points 
and  newly  constructed  mid-point  to  predict  values  of  the  Brownian  motion  at  the  one-fourth 
and  three- fourths  points  of  the  interval.  Random  values,  with  appropriate  error  variances, 
are  then  added  to  the  pi-edictious  at  each  of  these  new  points.  The  critical  observation  to 
be  made  here  is  that,  since  the  Brownian  motion  process  is  a  Markov  process,  its  value  at 
the  one- fourth  point,  given  the  values  at  the  initial  point  and  mid-point  is  independent  of 
the  process  values  beyond  the  mid-point,  in  particular  the  values  at  the  three- fourths  and 
end-points  of  the  interval.  Obviously,  it  is  also  the  case  that  the  value  at  the  three-fourths 
point  is  independent  of  the  values  at  the  initial  and  one-fourth  points,  given  the  values  at 
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the  mid-point  and  final  point.  TJie  consequence  of  this  observation  is  that  the  generation  of 
the  one-fourth  and  three-fourths  point  values  can  be  carried  out  in  a  parallel  and  decoupled 
fashion:  for  each  we  only  need  the  two  previously  generated  points  closest  to  it  (e.g.  the 
initial  and  mid-points  for  the  prediction  of  the  one- fourth  point)  plus  a  random  value  to 
be  added  to  this  prediction  that  is  independent  of  the  random  value  to  be  added  at  the 
other  point.  Thus,  each  of  these  two  steps  looks  identical  in  structure  to  the  original  step 
which  generated  the  mid-point  from  the  two  boundary  points.  In  addition,  we  see  that  the 
Markov  property  of  Brownian  motion  allows  us  to  iterate  this  process,  generating  values 
at  increasingly  dense  sets  of  dyadic  points  in  the  interval.  At  each  level  in  this  procedure 
we  generate  values  at  the  mid-points  of  all  neighboring  pairs  of  points,  in  parallel  and 
independently  of  previously  generated  points.  The  structure  of  all  of  these  steps  is  identical 
to  that  used  to  generate  the  original  mid-point. 

There  are  several  important  observations  to  be  made  about  the  preceding  development. 
The  first  is  that,  by  linearly  interpolating  at  each  level  in  this  procedure,  as  illustrated  in 
Figure  1,  a  sequence  of  continuous  approximations  of  a  Brownian  motion  is  constructed,  and 
the  statistics  of  these  approximations  converge  to  those  of  a  Brownian  motion  [27].  Indeed, 
this  sequence  of  linear  spline  approximations  can  be  interpreted  exactly  as  a  non-orthogonal 
multiscale  approximation  using  as  the  scaling  function  the  triangular  “hat”  function  [57] 
which  is  the  integral  of  the  Ilaar  wavelet  [-32].  Second,  as  we  will  see,  the  structure  of 
this  mid-point  deflection  construction  fits  precisely  into  the  multiscale  modeling  framework 
developed  in  [16,  17,  18],  and  corresponds  simply  to  a  particular  choice  of  the  parameters 
in  the  multiscale  model.  Moreover,  this  concept  generalizes,  allowing  us  to  show  that 
all  reciprocal  and  Markov  processes  in  one  dimension*  can  be  represented  by  multiscale 
stochastic  models  in  a  similar  way.  Thus,  in  one  dimension  we  will  show  that  the  class  of 
processes  realizable  via  multiscale,  scale-recursive  models  is  at  least  as  rich  as  the  class  of 
all  Markov  and  reciprocal  processes.  In  fact,  as  we  will  illustrate,  it  is  significantly  richer 
than  this. 

Furthermore,  and  most  significantly,  these  same  ideas  can  be  extended  to  multidimen¬ 
sional  processes.  In  particular,  we  show  how  a  generalization  of  the  mid-point  deflection 
concept  to  a  “mid-line”  deflection  construction  can  be  used  to  represent  all  2-D  Markov 
random  fields  with  multiscale  models.  In  particular,  the  key  to  our  multiscale  represen¬ 
tations  in  one  or  two  dimensions  is  a  partitioning  of  the  domain  over  w'hich  the  process 
is  defined  so  that  the  coarse-to-fine  construction  of  the  process  can  proceed  independently 
in  each  subdomain.  Markovianity  plus  knowledge  of  the  process  on  the  boundaries  of  the 
subdomain  partition  make  this  possible.  The  fundamental  difference,  however,  between  the 
1-D  and  2-D  cases  is  due  to  the  fact  that  boundaries  in  correspond  to  curves  or  in  2^ 
to  se/.sof  connected  lattice  sites,  as  opposed  to  pairs  of  points  in  one  dimension.  Because  of 
this  difference,  exact  multiscale  representations  of  MRF's  defined  over  a  subset  of  2'^  have 
a.  dimension  which  varies  from  scale  to  scale,  and  wdiich  depends  on  the  size  of  the  domain 
over  which  the  MRF  is  defined. 

As  a  consequence,  in  addition  to  the  exact  representations,  we  will  introduce  a  family  of 

'‘Note  that  this  also  includes  all  so-called  higher  order  Markov  and  reciprocal  processes.  For  example, 
a  second  order  Markov  process,  i.e.  one  for  which  the  value  of  the  process  i(f)  at  time  t  depends  on  both 
z{t  —  1)  and  z{t  —  2),  can  be  represented  as  a  recto?- first  order  Markov  process  [j(t),;(t  —  1)]*^,  which  can 
then  be  represented  in  the  manner  developed  in  Section  3. 
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approximate  representations  for  Gaussian  MRF’s  (GMRF’s)  based  on  wavelet  transforms. 
As  we  have  indicated,  maintaining  complete  knowledge  of  a  process  on  2-D  boundaries  leads 
to  models  of  scale-varying  dimension,  which  can  become  prohibitively  large  for  domains  of 
substantial  size.  On  the  other  hand,  at  coarser  scales,  it  would  seem  reasonable  to  keep 
only  coarse  approximations  to  these  boundary  values,  and  this  leads  naturally  to  the  use 
of  a  multiscale  change  of  basis  for  the  representatioji  of  the  values  of  a  2-D  process  along 
each  1-D  boundary.  That  is,  through  our  mid-line  deflection  based  models,  we  are  led  to 
the  idea  of  using  one- dimensional  wavelet  transforms  in  the  representation  of  the  values  of 
a  two-dimensional  GMRF. 

If  wavelet  coefficients  at  all  scales  are  kept  along  every  subdomain  boundary  for  the 
GMRF,  we  simply  have  an  exact  representation  for  the  GMRF  in  a  transformed  basis. 
However,  from  wavelet  analysis  of  1-D  stochastic  processes,  such  as  in  [20,  6,  61,  32],  we 
know  that  the  wavelet  transform  can  achieve  high  levels  of  scale- to-scale  decorrelation. 
This  result  has  led,  in  one  dimension,  to  approximate  models  and  algorithms  [20,  66]  which 
neglect  the  residual  correlation,  and  it  suggests  here  the  idea  of  keei)ing  only  a  subset  of  the 
wavelet  coefficients  representing  each  boundary.  While  this  approach  would  yield  a  coarse 
approximation  of  the  GMRF  along  each  boundary,  as  we  move  to  finer  scales,  the  length 
of  each  boundary  is  successively  halved,  so  that  “coarse”  approximations  at  fine  scales  are, 
in  fact,  increasingly  fine  (and,  indeed,  at  fine  enough  scales  represent  complete  wavelet 
descriptions).  The  result  is  a  family  of  models,  ranging  from  those  which  keep  only  the 
coarsest  wavelet  coefficients  along  each  1-D  boundary  to  the  exact  model  which  keeps  them 
all.  This  family  of  approximate  representations  allows  one  to  tradeoff  the  complexity  anil 
accuracy  of  the  representations,  while  also  providing  a  framework  for  the  development  of 
extremely  efficient  estimation  and  likelihood  calculation  algorithms. 

We  demonstrate  our  framework  for  wavelet-based  approximate  representation  of  Gaus¬ 
sian  MRF’s  in  the  context  of  natural  texture  I’epresentation.  In  particular,  classes  of 
GMRF's  have  been  widely  used  to  represent  natural  textures  in  the  context  of  segmen¬ 
tation  and  anomaly  detection  applications  [13,  22,  23,  24,  36,  47,  48],  and  we  illustrate  how 
these  models  can  be  approximated  in  our  multiscale  framework.  In  addition,  we  illustrate 
how  the  fidelity  of  the  approximation  varies  with  the  characteristics  of  the  GMRF  being 
approximated  and  with  the  complexity  of  the  approximate  representation. 

This  paper  is  organized  as  follows.  Section  2  describes  the  class  of  multiscale  stochastic 
models  that  we  use.  Section  3  develops  the  details  of  the  representation  of  Brownian  motion 
discussed  above,  and  generalizes  this  idea  to  allow  the  representation  of  all  1-D  Markov  and 
reciprocal  processes.  Section  4  then  describes  how  these  ideas  can  be  further  generalized 
to  provide  exact  and  approximate  representations  of  MRF’s.  Section  5  illustrates  how  the 
approximate  models  can  be  used  to  represent  GMRF  texture  models.  In  our  opinion,  one 
of  the  conclusions  that  can  be  drawn  from  these  results  is  that  this  multiscale  modeling 
framework  holds  substantial  promise  as  an  alternative  to  MRF’s  as  it  possesses  advantages 
both  in  efficient  optimal  algorithms  and  in  the  expressive  power  of  the  framework.  A 
number  of  interesting  and  substantive  problems  remain  to  be  investigated  —  such  as  the 
use  of  alternatives  to  wavelet  expansions  such  as  wavelet  packets  [62]  or  lapped  orthogonal 
transforms  [4.5,  46]  —  and  several  of  these  problems,  as  well  as  the  conclusions  of  this  paper, 
are  presented  in  Section  6. 
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2  Multiscale  Stochastic  Models 


In  this  section  we  describe  the  classes  of  multiscale  stochastic  models  to  be  used  in  this 
paper.  A  class  of  models  for  Gaussian  processes  is  described  first,  followed  by  a  simple 
generalization  allowing  for  more  general  (non-Gaussian)  processes.  For  simplicity,  in  this 
section  we  introduce  the  basic  structure  and  form  of  our  models  in  the  context  of  rej)- 
resenting  1-D  signals  and  processes.  The  extension  of  the  models  to  2-D  is  conceptually 
straightforward,  adding  only  notational  and  graphical  complexity,  and  thus  we  defer  the 
introduction  of  this  extension  until  Section  4,  where  it  is  needed. 

2.1  Gaussian  Multiscale  Models 

The  models  presented  in  this  section  were  recently  introduced  in  [16, 17, 18, 19]  and  describe 
multiscale  Gaussian  stochastic  processes  indexed  by  nodes  on  the  dyadic  tree®  shown  in 
Figure  2.  Different  levels  of  the  tree  correspond  to  different  scales  of  the  process.  In 
particular,  the  2”’^“^  values  at  the  m*’''  level  of  the  tree  are  interpreted  as  information  about 
the  scale  of  the  process,  where  the  notion  of  “information”  at  this  point  is  abstract.  For 
instance,  values  of  the  process  at  level  m  may  correspond  to  averages  of  pairs  of  values  at 
level  m  +  1.  In  this  case,  one  could  interpret  the  values  of  the  multiscale  process  as  scaling 
coefficients  in  a  Haar  wavelet  representation  of  the  process  at  the  finest  scale.  However, 
there  are  many  other  possible  interpretations  of  the  information  represented  at  each  level 
in  the  tree.  For  example,  values  of  the  multiscale  process  at  a  certain  level  could  also 
correspond  to  new  details  of  the  process  not  present  at  coarser  resolutions.  In  this  case,  the 
process  values  would  be  interpreted  as  the  wavelet  coefficients  in  a  wavelet  representation 
of  a  1-D  function  or  sequence.  Alternatively,  the  values  of  the  process  at  different  levels 
may  correspond  to  decimated  versions  of  the  process  at  the  finest  scale.  In  addition,  we 
may  wish  to  maintain  a  set  of  values  at  each  node  of  the  tree,  i.e.  a  “state”  vector,  in  order 
to  capture  the  scale-to-scale  memory  in  a  process.  For  example,  one  can  construct  higher 
order  models  in  which  the  state  includes  both  scaling  and  wavelet  coefficients  or  sets  of 
values  of  decimated  versions  of  the  process.  As  we  will  see,  our  multiscale  representations 
of  reciprocal  processes  and  MRF’s  have  a  natural  interpretation  in  terms  of  decimated 
versions  of  the  processes  being  modeled,  although  they  can  also  be  interpreted  in  terms  of 
scaling  coefficients  corresponding  to  particular  non-orthogonal  multiscale  representations. 
Moreover,  these  models  will,  in  fact,  require  the  use  of  state  vectors  at  each  node  on  the 
tree. 

We  denote  nodes  on  the  tree  with  an  abstract  index  s,  and  define  an  upward  shift 
operator  7  such  that  57  is  the  parent  of  node  s,  as  illustrated  in  Figure  2.  Note  that  the 
operator  7  is  not  one-to-one,  but  is,  in  fact,  two-to-one  since  each  parent  has  two  offspring. 
Also,  we  define  the  scale  of  node  s,  i.e.  the  level  of  the  node,  as  ni(s).  The  stochastic  tree 
process  x(s)  £  is  then  described  via  the  following  scale- recursive  dynamic  model  driven 
by  white  noise: 


X  ( s )  =  A(  .s  )x(s'))  +  B{s)w(s)  ( 1 ) 

^As  we  will  see  in  Section  4  and  as  is  illustrated  in  Figure  12,  the  extension  to  2-D  signals  involves,  in 
essence,  replacing  the  dyadic  tree  by  a  quadtree. 
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m=  1 


m  =  2 

m  =  3 

m  =  4 


Figure  2:  The  state  vectors  in  multiscale  stochastic  models  are  indexed  by  the  nodes  of  a 
dyadic  tree.  The  tree  is  a  set  of  connected  nodes,  in  which  each  node  has  two  offspring. 
The  parent  of  node  s  is  denoted  and  the  scale,  or  level,  of  node  s  is  denoted  by  in(s). 
Multiscale  stochastic  processes  defined  on  higher-order  trees  (see  Figures  7  and  12)  will  also 
be  used. 

under  the  following  assumptions®: 


.ro  ~  A7(0,Fo)  (2) 

w(s)  ~  A^(0,/)  (3) 

where  u){.s)  G  and  A{s)  and  B{s)  are  matrices  of  appropriate  size.  The  state  variable 
Xo  at  the  root  node  of  the  tree  provides  an  initial  condition  for  the  recursion.  The  driving 
noise  tcis)  is  white,  i.e.  to(.s)  and  ^re  nncorrelated  if  -s  a.  and  is  nnrorrebvted 

with  the  initial  condition.  Interpreting  each  level  as  a  representation  of  one  scale  of  the 
process,  we  see  that  (1)  describes  the  evolution  of  the  process  from  coarse  to  fine  scales. 
The  term  yl(.s).x-(.s7)  represents  interpolation  or  prediction  down  to  the  next  level,  and 
B(.9)w{s)  represents  new  information  added  as  the  process  evolves  from  one  scale  to  the 
next.  The  choice  of  the  parameters  .T(s)  and  B{.s)  and  their  dependence  (if  any)  on  the 
node  s,  depends  upon  the  particular  application  and  process  being  modeled.  For  example, 

®Tlie  notation  x  ~  Af{m,  P)  means  the  random  vector  x  is  normally  distributed  v^ith  mean  rn  and  variance 
P.  We  also  sometimes  use  the  notation  px(X)  =  A^(A';  m,  P)  to  denote  the  same  thing. 
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if  a'(5)  is  scalar  and  yl(5)  =  1,  the  values  of  x{.s)  at  any  particular  scale  have  a  natural 
interpretation  as  providing  a  piecewise-constant  (i.e.  Haar)  approximation  to  a  signal,  and 
the  B(.s)w{s)  terms  then  represent  additional  detail  added  in  progressing  from  coarse  to  flue 
scales  [10,  17,  18,  19].  If,  in  addition,  the  magnitude  of  this  detail  decreases  geometrically 
with  scale  —  i.e.  if  B(s)  =  where  /j.  is  a  scalar  constant  —  then  the  process  so 

represented  has  fractal  characteristics  [17,  18,  20].  In  fact,  in  general,  dependence  of  /l(.s) 
and  B{s)  on  scale  m{s)  allows  one  to  model  self-similar  statistical  structures  or  the  pi’esence 
of  dominant  scales  [20],  while  allowing  general  dependency  on  the  node  s  provides  additional 
flexibility  including  the  ability  to  capture  non-stationarities,  localized  events,  etc.  In  the 
context  of  this  paper,  as  we  will  see,  the  parameters  of  the  model  (1)  are  determined  in  a 
constructive  fashion  in  order  to  represent  the  reciprocal  process  or  MRF  of  interest. 

An  extremely  important  property  of  the  scale- recursive  model  (1)  is  that  not  only  is 
it  Markov  from  scale- to-scale  but  it  also  is  Markov  with  respect  to  the  partial  ordering 
defined  by  the  dyadic  tree  structure.  That  is,  conditioned  on  the  vabte  at  a  parent  node, 
the  evolution  of  the  process  in  the  subtree  descendant  from  that  node  is  independent  of 
everything  outside  that  subtree.  This  fact  implies  that  there  are  extremely  efficient  and 
highly  parallelizable  algorithms  for  optimal  estimation  and  other  signal  processing  tasks 
based  on  this  model. 

In  particular,  as  developed  in  [2,  16,  17,  18,  19],  there  is  an  extremely  efficient  algorithm 
for  optimal  estimation  based  on  noisy  measurements  y{s)  6  of  the  process  of  the  form: 

y{s)  =  Cfs^rfs)  -b  u(s)  (4) 

where  v(s)  ~  yV(0,i?(s))  and  the  matrix  C7(s)  can  sjrecify,  in  a  very  general  way,  measure¬ 
ments  taken  at  diflerent  times  or  spatial  locations  and  at  different  scales.  The  structure  of 
the  algorithm  for  this  problem  consists  of  two  scale-recursive  sweeps  (a  flne-to-coarse  sweep 
followed  by  a  coarse-to-fine  recursion),  each  of  which  follows  the  structure  of  the  dyadic  tree 
so  that  there  is  substantial  parallelism  and  efficiency.  For  example,  the  extension  of  this  al¬ 
gorithm  to  2-D  and  quadtrees  is  applied  in  [42]  to  develop  a  new  scale-recursive  approach  to 
dense  motion-field  estimation  in  image  sequences  that  is  considerably  faster  than  previously 
developed  algorithms.  In  addition,  a  related  algorithm,  which  is  also  highly  parallelizable, 
can  be  used  to  compute  the  likelihood  of  the  measurements  (log  of  the  conditional  proba¬ 
bility  density)  given  the  model  parameters  [4.3].  This  algorithm  can  be  used,  together  with 
the  results  presented  here,  for  texture  identification  and  segmentation,  for  example.  An  im¬ 
portant  point  about  these  algorithms,  which  is  of  particular  significance  for  2-D  processing, 
is  that  they  are  recursive  and  not  iterative,  and  in  fact  have  constant  complexity  per  data 
point  or  pixel.  This  is  in  sharp  contrast  to  the  usual  iterative  algorithms  associated  with 
the  processing  of  MRF’s  [3.3j. 

2.2  General  Multiscale  Models 

As  we  indicated  in  the  preceding  section,  a  basic  property  of  the  model  (1)  -  (3)  is  the 
Markovianity  of  the  state  with  respect  to  the  ordering  structure  defined  by  the  dyadic  tree. 
More  precisely,  two  states  are  conditionally  independent,  given  a  state  on  the  path  between 
them: 

fA:(ai),a,-(s2)b’(«3),S3ersj,s2  1  >  "^3  G  — 
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7\-(si)|i-(s3),a3€r,j,«2(^''»ll'^*S3’'-^3  €  )Px(s2)|.r(a3),S36r,j,aj  G  I'si  ,52  ) 

where  is  defined  as  the  set  of  nodes  in  the  shortest  i>ath  between  two  nodes  s  and  <t  (but 
not  including  s  or  a).  By  starting  with  this  assumption,  and  then  allowing  the  probability 
density  function  (pdf)  for  the  state  at  a  particidar  node  given  its  parent  and  the  pdf  for  the 
state  at  the  root  node  to  be  arbitrary,  we  obtain  a  much  wider  class  of  processes  than  that 
given  by  (1)  -  (3).  Indeed,  we  define  a  Markov  tree  process  as  any  set  of  random  variables 
which  are  indexed  by  the  nodes  of  a  tree  and  have  a  joint  pdf  satisfying  (5). 

Markov  tree  processes  are  naturally  defined  by  specifying  the  parent-offspring  condi¬ 
tional  pdf’s,  along  with  a  pdf  for  the  state  at  the  root  node  of  the  tree.  A  simple  example 
of  a  stochastic  process  in  this  more  general  class  is  the  following  discrete-state  stochastic 
process  x(s)  G  {(),  1,  •  •  • ,  A)  with  parent-offspring  conditional  probability  mass  functions 
given  by: 


2A’(s)|a;(s7)  (-^  s|  A 

57) 

7>aTo(A’o) 


J  ^m(s)  if  A3  —  A  3-), 

\  (l-^m(s))/i  ifA'3  5^A'3^ 

1/(1 +  1),  for  AoG  {0,l,---,i} 


(6) 

(7) 


where  **  number  between  0  and  1  which  may  vary  with  scale  m(s).  A  class  of 

processes  with  this  structure  and  defined  on  a  quadtree  has  been  proposed  by  Bouman 
for  segmentation  applications  [7,  8,  9].  We  will  see  several  other  examples  of  such  general 
Markov  tree  processes  in  the  next  section. 

The  property  (.5)  implies  that  the  tree  processes  are  Markov  iji  scale,  from  coarse-to-fiue. 
In  particular,  we  denote  the  set  of  states  at  a  given  scale,  say  mo,  as  =  {;r(.s)|n?(.«)  = 
mo).  Then,  using  (.5),  the  states  at  a  particular  scale,  given  the  states  at  the  previous  level, 
are  independent  of  all  coarser  scales: 

<  ??r2)  =  7V”2|3:>»2-i(A'”"^|A'”'2-^)  (8) 

Indeed,  (.5)  implies  the  conditional  pdf  of  the  state  at  node  s,  given  the  states  at  all  previous 
scales,  depends  only  on  the  state  at  the  parent  node  57: 


Pa:(s)b’(cr),m(iT)<m(s)(A3|ArCT,  mfcr)  <  77r(s))  Pa^(5)|i’(s7)(A*s|A'3..^. )  (9) 

In  addition,  (.5)  implies  that  such  a  process  can  also  be  viewed  as  a  Markov  random  field 
on  the  tree.  In  particular,  define  the  neighborhood  set  Ds  of  node  s  as  its  three  nearest 
neighbors  on  the  tree'^  (i.e.  the  parent  and  two  offspring  nodes).  The  general  multiscale 
process  described  above  is  then  a  Markov  random  field  on  the  tree  in  the  sense  that: 

Fa^(s)|a'(<T),cr^s(  As|A(7,  (7  ^  .s)  =  ib-(s)|a'((T),<TgDs(  A*s|  A^,  O'  G  Z?3  )  (10) 

The  most  general  class  of  joint  probability  distribution  functions  which  lead  to  conditional 
distributions  satisfying  (10)  is  given  by  the  Hammersley-Clifford  theorem  [4].  Our  focus 
here  is  on  processes  which  also  satisfy  the  one-sided  Markov  property  (9).  because  this  class 
of  processes  leads  naturally  to  efficient  scale-recursive  algorithms. 

'Obvious  modifications  of  the  neighborhood  set  must  be  made  for  the  root  node  at  the  top  of  the  tree, 
which  has  no  parent,  and  the  nodes  at  the  finest  level  of  the  tree,  which  have  no  offspring. 
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Finally,  we  stress  that  while  (5)  implies  tliat  a  tree  process  is  an  MRF  with  respect  to 
the  nearest  neighbor  set  on  the  tree,  the  set  of  states  x”*  at  scale  m,  viewed  as  a  sequence  of 
length  2"^"^  is  not  Markov  for  an  arbitrarily  chosen  set  of  parent-offspring  pdf’s.  This  point 
can  be  appreciated  by,  for  example,  computing  the  joint  pdf  for  the  four  values  at  the  third 
level  of  the  multiscale  process  given  by  (6),  and  directly  checking  the  conditions  required 
for  Markovianity  of  the  single  level  sequence®.  However,  as  we  show  in  the  next  section, 
the  parent-offspring  conditional  pdf’s  can  be  chosen  such  that  the  single  level  process  r”' 
is  Markov  in  the  usual  sense,  and  in  fact  such  that  the  finest  level  of  the  tree  process  can 
be  used  to  represent  any  1-D  Markov  or  reciprocal  process,  with  higher  levels  in  the  tree 
corresponding  to  representations  of  the  process  at  coarser  resolutions. 

3  Representation  of  1-D  Reciprocal  and  Markov  Processes 

In  this  section  we  describe  the  basic  properties  of  reciprocal  processes  in  one  dimension,  in¬ 
troduce  and  develop  representations  of  reciprocal  processes  in  terms  of  multiscale  stochastic 
models,  and  present  several  examples. 

3.1  1-D  Reciprocal  Processes 

A  reciprocal  process  is  nothing  more  than  a  first-order  MHF  on  the  real  line.  More  formally, 
a  stochastic  process  z{t),t  £  TZ  is  said  to  be  reciprocal®  (or  bilateral  Markov,  two-sided 
Markov  or  non-causal  Markov)  if  it  has  the  property  that  the  probability  distribution  of 
a  state  in  any  open  interval  (Ti,T2),  conditioned  on  the  boundary  states  z(2\),z{T2)  is 
independent  of  states  outside  of  the  interval  [29,  38].  That  is,  for  t  €  (Ti,T2): 

Pz{t)\z(T),Te(Ti,T2)'--(Zt\ZT,T  e  iTi,T2f)  = 

Pz(t)\z(Ti)MT2)i^1:\ZTi-,  ZT2)  (11) 

where  (T'i,T2)'^  denotes  the  complement  of  the  open  interval  (T1/T2).  Reciprocal  processes 
defined  on  the  integers  Z  satisfy  the  same  property  with  the  continuous  interval  {2\,T2) 
replaced  by  the  discrete  interval  {Ti  -f-  1,  Ti  -f  2,  •  •  ■ ,  r2  —  1}: 

Pz{t)\z(T),T^{Tl-^\,-,T2-l-YiZt\ZT,  T  6  {Tf  -h  1,  ■  •  •  ,  r2  -  1}'^)  = 

Pz(t)\z(Ti),z(T2)i^t\^Tx,  ZT2)  (12) 

Reciprocal  processes  are  closely  related  to  the  class  of  Markov  processes.  A  process 
;;(f)  on  TZ  or  Z  is  Markov  if  past  and  future  values  of  the  state  are  independent  given  the 
present.  This  means  that  for  #2  <  fs: 

Pz(tz)\z(ti),tl<t2^^i3\Ztx-,h  <h)  =  P-(t:^)\;(t2)iZti\Zt2)  (13) 

As  shown  in  [l],  if  a  process  is  Markov  then  it  is  also  reciprocal.  On  the  other  hand, 
reciprocal  processes  are  not  necessarily  Markov  [29],  although  one  can  show  that  essentially 
all  stationary  Gaussian  reciprocal  processes  are  Markov  [40]. 

®The  process  is  Markov  only  if  0m(s)  =  l/(-t  +  !)•  In  this  case,  the  values  of  the  process  at  any  level  in 
the  tree  are  independent  of  one  another. 

®The  discussion  here  refers  only  to  first-order  reciprocal  processes.  Extension  to  higher-order  processes 
is  straightforward  [29]. 
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3.2  Exact  Multiscale  Representations  of  1-D  Reciprocal  Processes 

In  the  introduction  we  described  a  construction  of  a  Brownian  motion  b(t)  over  an  interval, 
say  [0, 1],  via  mid-point  deflection  [27].  As  we  noted,  this  corresponds  precisely  to  one  of 
the  Gaussian  multiscale  stochastic  models  described  in  .Section  2.  To  see  this,  consider  the 
following  process.  At  the  coarsest  level,  the  initial  state  .^•o  is  a  three-dimensional  vector 
whose  pdf  is  given  by  the  joint  pdf  for  the  values  of  a  Brownian  motion  at  the  initial,  middle 
and  final  points  of  the  interval: 

6(0) 

xo  =  6(0.5)  ~A^(0,Po) 

.  ^(1)  . 

■  0  0  o' 

Po  =  0  0.5  0.5 

0  0.5  1 

where  we  have  used  the  facts  that  6(0)  =  0,  b{t)  is  an  independent  increments  process,  and 
for  ii  <  t2,  b(t2)  -  b(ti)  ~  A/'(0,t2  -  h)- 

Choosing  a  value  for  .ro  as  a  sample  from  the  distribution  (14)  corresponds  to  the  first 
two  steps  in  the  mid-point  deflection  construction  of  Browniaji  motion,  the  first  step  being  a 
choice  for  the  two  end-point  values  6(0),  6(1),  and  the  second  step  being  construction  of  the 
mid-point  6(0.5).  The  next  step  in  the  mid-point  deflection  construction  is  the  specification 
of  values  for  the  Brownian  motion  at  the  one-fourth  and  three- fourths  points.  In  the  context 
of  our  multiscale  modeling  framework,  we  define  two  state  vectors  at  the  second  level  of 
the  dyadic  tree  in  Figure  2,  each  again  a  3-tuple.  The  state  on  the  left  represents  the 
values  of  the  Brownian  motion  at  the  initial,  one-fourth  and  middle  points  of  the  interval, 
[6(0), 6(0.25), 6(0.5)],  and  the  state  on  the  right  represents  the  corresponding  values  in  the 
right  half-interval,  [6(0.5),  6(0.75),  6(1)].  The  sample  at  the  quarter  point  is  given  by  linear 
interpolation  of  6(0)  and  6(0.5),  plus  a  Gaussian  randojn  variable  with  variance  equal  to  the 
variance  of  the  eri'or  in  this  prediction: 

6(0.2.5)  =  i(6(0)  + 6(0.5)) -)-e(0.25)  (16) 

e(0.25)  ~  A^(0, 0.125)  (17) 

Likewise,  6(0.75)  is  chosen  by  averaging  the  end  points  of  the  right  half-interval,  6(0.5)  and 
6(1),  and  adding  in  a  random  value,  independent  of  the  deflection  term  used  to  create  the 
sample  at  the  one-fourth  point: 

6(0.75)  =  i(6(0.5)  + 6(1)) +  6(0.75)  (18) 

e(0.75)  ~  A^(0, 0.125)  (19) 

The  above  construction  of  6(0.25)  and  6(0.75)  is  precisely  the  same  as  the  mid-point  de¬ 
flection  construction  of  these  values.  Values  of  the  process  at  successively  finer  sets  of 
dyadic  points  are  generated  in  the  same  way.  At  the  scale,  the  values  of  the  process 
at  t  =  k/2"^,k  =  0, 1,  •  •  •,2’”  are  represented  with  2”^~^  state  vectors,  each  containing  the 


(14) 

(15) 
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b(t) 


Scale 


m=  1 


m  =  2 
m  =  3 


Figure  3:  The  state  vectors  for  the  first  three  levels  of  a  multiscale  model  representing  Brow¬ 
nian  motion,  6(t),  are  illustrated.  At  the  first  level,  the  state  is  the  vector  [h(0),  6(0.5),  h(  1)], 
which  is  indicated  by  the  three  points  at  m  =  1  surrounded  by  an  ellipse.  The  points  are 
placed  directly  below  the  points  t  —  0, 0.5  and  t  =  1  on  the  graph  above  to  indicate  that 
the  state  of  the  multiscale  process  at  the  first  level  consists  of  the  values  of  the  Brownian 
motion  at  those  three  points.  Likewise,  at  lower  levels,  the  states  are  indicated  by  sets  of 
three  points  surrounded  by  ellipses,  with  the  horizontal  location  of  the  points  in  correspon¬ 
dence  with  time  indices  in  the  the  graph  at  the  top.  At  the  level,  there  are  state 
vectors,  each  of  which  consists  of  the  values  of  b(t)  at  three  consecutive  dyadic  points,  and 
which  together  represent  the  values  of  the  Brownian  motion  at  2”*  +  1  distinct  points  on 
the  interval  [0, 1].  The  multiscale  representation  for  Brownian  motion  can  be  generalized  to 
the  class  of  1-D  reciprocal  process,  which  contains  the  class  of  1-D  Markov  processes, 

values  of  the  process  at  three  points,  as  shown  in  Figure  3.  At  any  level,  each  state  is  a 
linear  function  of  its  parent,  plus  an  independent  noise  term.  Thus,  this  construction  fits 
precisely  into  the  multiscale  modehng  framew'ork  given  by  (1)  -  (3)  (see  Section  3.3  for  the 
precise  formulae  for  A(i!)  and  B(s)). 

Representation  of  more  general  1-D  reciprocal  processes  via.  multiscale  models  is  a  simple 
extension  of  the  above  idea.  To  construct  a  multiscale  model  for  a  particular  reciprocal 
process  z(t),  t  6  [0, 1],  start  by  choosing  the  state  at  the  coarsest  level  as  a  sample  from  the 
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joint  distribution: 


Pz{0),z(Q.^),z(l){^0i  Zo.5,  Zi).  (20) 

This  generalizes  the  choice  in  ( 14)  in  wliicli  the  state  at  the  top  level  is  chosen  using  the  Gaus¬ 
sian  distribution  corresponding  to  a  Brownian  motion.  The  two  state  vectors  at  the  second 
level  are  again  the  three-dimensional  vectors  [r(0),  ^(0.25),  ^(0.6)]  and  [r(0..5),  r(0.75 ),  z(  1 )], 
where  values  for  the  half-interval  mid-points  are  chosen  as  samples  from  the  conditional  dis¬ 
tributions: 

Pz(0.25)\z(0),z{0.5)(Zo.25\Zo,  Zq.s)  (21) 

14(0.75)1^(0.5), s(1)(-^0.7.5|-^0.5j  Zi  )  (22) 

Since  the  process  is  reciprocal,  z(0.25)  and  a:(0.75)  are  conditionally  independent  given  the 
state  at  the  first  level,  and  thus  the  modeling  structure  fits  precisely  into  the  more  general 
non-linear  model  class  described  in  Section  2.2. 

The  construction  above  assumes  that  the  process  is  defined  over  a  confinuous  inter¬ 
val.  In  practice,  we  are  typically  concerned  with  processes  z{t)  on  a  discrete  interval, 
t  e  {0, 1, •  •  •  ,T'}.  If  r  =  2^  for  some  integer  iV,  then  we  can  use  essentially  the  same 

construction  as  for  the  continuous  case  above.  Specifically,  xq  =  [z(0),  z(T/2},  z{T)]  is  a 

random  vector  chosen  from  the  appropriate  distribution  for  the  process  of  interest.  The 
states  at  the  second  level  are  lz(0),z(T/4),z(T/2)]  and  [z(T/2),z(3T/4),z(T}],  with  the 
half-interval  mid-points  again  chosen  using  the  appropriate  distribution.  Since  there  are 

only  a  finite  number  of  points  in  the  discrete  proce,ss,  only  a  finite  number  of  levels  are 

needed  to  exactly  represent  it.  In  particular,  with  T  =  2^^,  i\^  levels  are  required. 

There  are  several  observations  to  be  made  about  the  continuous  and  discrete- time  con¬ 
struction  we  have  just  described.  The  first  is  that  there  is  no  fundamental  difficulty  in 
choosing  a  point  other  than  the  mid-point  at  ea.ch  level  in  these  constructions.  For  exam¬ 
ple,  in  the  construction  of  Brownian  motion,  starting  from  the  initial  set  of  points  given  in 
(14),  we  could  next  generate  any  pair  of  points  on  either  side  of  0.5,  e.g.  6(0.1)  and  6(0.7). 
However,  the  regular  structure  implied  by  the  choice  of  mid-points  may  be  of  some  value 
for  processes  such  as  Brownian  motion  which  have  stationary  increments,  as  they  lead  to 
models  in  which  the  model  parameters,  such  /l(.s)  and  B(s)  in  (1)  -  (.3)  have  very  simple 
and  regular  characterizations  as  a  function  of  node  s  and  scale  7?i(s)  (see,  for  example, 
(46),(47)).  This  regularity  in  turn  leads  to  simplifications  in  the  structure  of  algorithms 
for  estimation  and  signal  processing,  requiring  fewer  distinct  gains  to  be  calculated  and, 
if  parallel  implementation  is  considered,  allowing  SIMD  (single  instruction,  multiple  data) 
rather  than  MIMD  (multiple  instruction,  multiple  data)  implementations. 

Secondly,  in  discrete-time,  there  will  always  be  at  least  some  degree  of  irregularity  in 
the  mnltiscale  model  if  the  process  is  defined  over  t  E  {0. 1,  •  •  • .  T}  and  T  is  not  a  power  of 
two.  In  particular,  in  such  a  case  the  structure  of  the  tree  and/or  the  state  needed  in  the 
multiscale  representation  of  this  process  will  need  to  be  modified.  For  example,  consider 
a  process  defined  over  /  G  {0, 1,  •  •  • ,  10}.  In  this  case,  we  can  develop  a  model  of  the  type 
we  have  described  in  which  the  tree  is  of  non-uniform  depth  and  in  wdiich  we  do  not  have 
mid-point  deflection  at  some  nodes,  as  indicated  in  Figure  4a  (e.g.  in  the  generation  of  the 
value  at  t  =  3  given  values  at  0  and  5).  Alternatively,  as  shown  in  Figure  4b,  we  may  be 
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able  to  achieve  some  level  of  (and  perhaps  complete)  symmetry  by  genei-ating  more  than 
one  new  point  at  some  nodes  (e.g.  in  Figure  4b  we  generate  values  at  both  /  =  2  and  t  =  3 
given  values  at  0  and  5).  Obviously,  as  in  standard  discrete  signal  processing  applications 
in  which  the  FFT  is  to  be  used,  there  are  considerable  efficiencies  to  be  had  if  T  is  a  power 
of  2. 

Furthermore,  as  we  have  indicated  previously,  while  our  development  has  focused  on 
first-order  reciprocal  processes,  the  extension  to  higher-order  models  is  straightforward. 
Indeed,  a  A'^^-order  model  defined  on  t  6  {1,2,-  ••,  A'lF  +  1)},  where  T  is  a  power  of  2, 
can  be  accommodated  by  grouping  states  at  adjacent  points  into  sets  of  size  K.  The  states 
at  different  levels  of  the  tree  might  be  as  depicted  in  Figure  5  for  a  second-order  reciprocal 
process  defined  over  t  £  {0, 1,-  •  -  ,  18}.  Higher-order  models  can  equivalently  be  represented 
by  simply  redefining  the  state  of  the  process  z{t)  to  be  a  vector  of  appropriate  dimension. 

The  representations  we  have  introduced  to  this  point  have  obvious  and  substanf  ial  levels 
of  redundancy.  For  example,  the  value  of  z{T j'2)  appears  in  the  state  vector  at  both  nodes 
at  the  second  level  of  the  midtiscale  model  we  have  described  for  discrete-time  reciprocal 
processes.  More  generally,  at  the  level  of  the  model  for  such  a  process  there  are 
state  vectors  containing  a  total  of  3  X  2’"“^  values,  only  2™  -|- 1  of  which  are  distinct.  This 
redundancy  is  actually  of  minimal  consequence  for  estimation  algorithms  based  on  these 
models.  However,  it  is  also  easy  to  eliminate  the  redundancy  by  a  simple  modification  to  the 
construction  we  have  described.  In  particular,  we  may  generate  two  internal  points  between 
each  pair  of  points  at  each  stage  in  the  level-to-level  recursion,  yielding  a  four- dimensional 
state  vector.  For  example,  if  the  reciprocal  process  is  defined  over  t  £  {1, 2,  •  •  • ,  16),  then 
we  can  choose  the  non-redundaut  set  of  state  vectors  illustrated  in  Figure  6.  In  this  case, 
a  first-order  reciprocal  process  is  represented  by  a  process  with  a  four-dimensional  state, 
instead  of  the  process  with  a  three-dimensional  state  used  earlier.  In  general,  at  the 
level  of  such  a  representation,  there  are  2”*“^  state  vectors  representing  2™+^  distinct  values 
of  the  process.  Again,  in  the  situation  where  T  is  not  a  power  of  two,  some  irregularity  in 
the  structure  wall  be  introduced. 

Once  we  allow  ourselves  to  consider  such  variants  on  the  original  mid-point  deflection 
construction  in  which  more  than  one  new  point  is  generated  between  each  pair  of  previously 
constructed  points,  we  see  immediately  that  it  is  possible  to  generate  multiscale  representa¬ 
tions  on  trees  that  are  not  dyadic.  For  example,  a  q^^-order  tree  is  a  set  of  connected  nodes, 
such  that  each  node  has  q  offspring  as  in  Figure  7.  Consider  a  reciprocal  process  defined  on 
t  £  {0, 1,  •  •  • ,  3^}.  This  process  is  most  conveniently  represented  on  the  regular  structure  of 
a  third-order  tree,  as  shown  in  Figure  8.  This  flexibility  of  the  modeling  framework  allows 
the  possibility  of  considering  different  tradeoffs  in  terms  of  level  of  parallelization  and  com¬ 
putational  pow'er  of  individual  processors  when  implementing  estimation  algorithms  such 
as  those  in  [2,  16,  17,  18,  19,  42]. 

Finally,  it  is  of  interest  to  note  that  the  construction  w'e  have  described,  and  its  sev¬ 
eral  variants,  can  be  interpreted  as  a  non-iterative  Gibb’s  sampler.  The  Gibb’s  sampler 
introduced  in  [33]  is  an  iterative  algorithm  for  the  generation  of  sample  paths  of  Markov 
random  fields  on  a  discrete  lattice.  For  1-D  discrete-time  reciprocal  processes  this  proce¬ 
dure  reduces  to  using  the  nearest  neighbor  conditional  probability  functions  to  construct  a 
Markov  chain  which  has  an  asymptotic  distribution  equal  to  the  correct  joint  distribution  of 
the  process.  Specifically,  at  each  step  of  the  procedure  we  modify  the  current  sample  path 
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m=  1 


m  =  3 


Figure  4.  The  state  vectoi’s  are  shown  for  two  possible  luultiscale  representations  for  a 
reciprocal  process  which  is  defined  on  a  discrete  interval  of  the  form  {0. 1,  •  •  ■ ,  10}.  In  (4a), 
a  dyadic  tree  with  uniform  state  dimension,  but  non-uniform  depth  is  used,  whereas  in  (4b) 
a  dyadic  tree  of  uniform  depth  but  non-uniform  state  size  is  used. 
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Figure  5:  The  state  vectors  are  shown  for  a  multiscale  representation  of  a  second-order 
reciprocal  process  defined  on  a  discrete  interval.  In  this  case,  the  state  vectors  are  composed 
of  three  groups  of  two  points,  reflecting  the  fact  that  the  value  at  the  current  point,  say 
z{to)  is  indeirendelit  of  the  values  at  all  other  times,  given  the  pairs  of  nearest  neighbors 
to  “  1)5  ^(to  “  2)  and  z(to  -f  1),  z{to  -f-  2). 


Zk(i),  t  £  {0, 1,  •  •  -  ,7}  by  replacing  the  value  at  some  point  in  time,  say  to  with  a  random 
value  chosen  according  to  the  conditional  distribution  for  the  process  at  that  point  given 
the  current  values  of  the  sample  path  at  to  —  1  and  to  +  1.  By  cycling  repeatedly  through  all 
of  the  time  points,  the  sample  path  is  guaranteed  to  converge  to  one  with  the  correct  joint 
statistics.  The  procedure  is  conceptually  simple  but  computationally  intensive,  .since  the 
Markov  chain  requires  many  transitions  for  the  probability  function  to  converge.  In  con¬ 
trast,  in  our  construction,  we  successively  generate  samples  at  new  points  (e.g.  mid-points) 
conditioned  on  values  at  previously  generated  points,  which  are  not  nearest  neighbors  but 
rather  boundary  points  that  partition  the  time  interval  of  interest.  For  this  reason,  and 
since  we  begin  at  the  root  node  with  a  decimated  set  of  values  with  the  correct  distribution, 
we  are  guaranteed  that  at  each  stage  the  decimated  process  that  is  constructed  has  exactly 
the  correct  distribution.  Thus,  with  this  structure  we  visit  each  time  point  only  once  and 
construct  a  sample  path  non-iteratively. 

In  fairness,  an  important  point  to  note  here  is  that  if  a  reciprocal  process  is  specified 
directly  in  terms  of  a  Gibb’s  distribution  —  i.e.  in  terms  of  nearest  neighbor  energy  functions 
(see,  for  example,  [3.3])  —  then  the  calculation  of  the  nearest  neighbor  pdf’s  required  in  the 
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Figure  6:  The  state  vectors  are  shown  for  a  non-redundant  multiscaJe  representation  of 
a  1-D  reciprocal  process.  These  non-redundant  representations,  appropriately  generalized 
for  the  2-D  case,  are  useful  in  the  context  of  wavelet-based  approximate  representations  of 
Gaussian  MRF’s. 


q  offspring 


Figure  7:  The  g*^-order  tree  is  illustrated.  This  is  a  set  of  nodes,  each  of  which  has  q 
offspring. 
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Figure  8:  The  state  vectors  are  shown  for  a  uiultiscale  representation  on  a  third-order  tree. 

Gibb's  sampler  is  simple.  The  question  then  is  whether  it  is  also  simple  to  determine  the 
conditional  pdf’s  —  e.g.  the  pdf  for  z{T  12)  given  ^(0)  and  z{T)  —  needed  to  implement 
the  non-iterative,  multiscale  procedures  we  have  described.  In  general,  this  may  not  l)e  a 
straightforward  task,  since,  if  we  begin  with  a  Gibb’s  distribution  the  computation  of  such 
pdf’s  requires  explicit  calculation  of  quantities  related  to  the  so-called  partition  function 
[33],  which  can  be  quite  complex.  On  the  other  hand,  such  calculations  need  only  be 
done  once  to  construct  the  pdf’s  needed  for  the  multiscale  model.  In  addition,  as  we  have 
seen  for  Brownian  motion  and  as  we  now  illustrate  further,  in  many  cases,  including  all 
vector  Gauss-Markov  processes,  closed  form  expressions  can  be  derived  for  the  multiscale 
representations. 

3.3  Examples 

Ill  this  section  we  discuss  several  examples  of  reciprocal  processes  and  their  multiscale  rep¬ 
resentations.  The  first  examples  describe  multiresolution  models  for  general  vector  Gauss- 
Markov  processes  specified  in  state-space  form,  and,  in  particular  describe  this  construction 
in  detail  for  two  cases  corresponding  to  the  integral  of  white  noise  (i.e.  Brownian  motion) 
and  the  second  integral  of  white  noise.  These  examples  allow  us  to  illustrate  the  interpre¬ 
tation  of  these  multiresolution  models  as  providing  approximations  using  noii-orthogonal 
expansions.  In  particular,  our  model  for  Brownian  motion  corresponds  to  the  use  of  the  so- 
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called  “hat”  function  [57]  in  this  expansion,  leading  to  linear  interpolation  between  dyadic 
points,  while  the  model  for  the  integral  of  Brownian  motion  leads  to  a  multiresolution 
approximation  using  cubic  interpolation. 

The  second  part  of  this  section  presents  several  sijiiple  discrete-state  examples,  the  first 
of  which  investigates  general  A'^-state  Markov  chains  and  allows  us  to  make  contact  with 
the  models  used  in  [7,  8,  9]  for  segmentation  applications.  The  second  example  is  a  general 
two-state  process,  which  is  used  to  demonstrate  that  the  class  of  multiscale  models  is  in 
fact  far  richer  than  the  class  of  Markov  processes.  In  particular,  through  this  example 
we  gain  insight  into  the  very  particular-  conditions  that  the  parent-to-child  transition  pdf’s 
must  satisfy  in  order  for  the  finest  level  process  to  be  Markov.  This  analysis  suggests  that 
Markov  chains  are,  in  fact,  only  a  small  subset  of  the  processes  realizable  with  imdtiscale 
models  and,  in  particular,  directly  motivates  several  other  multiscale  mid-point  deflection 
processes  which  are  not  Markov. 

3.3.1  Gauss-Markov  Processes 

Consider  a  vector  Gauss-Markov  process  defined  on  the  interval  [0, 1]  and  given  by: 


z{t)  =  Ftz{t) -{■  Gti.i{t)  (23) 

where: 

^(0)  ~  .V(0,no)  (24) 

=  JS{t-T)  (25) 

E{fiit)z(0f}  =  0  (26) 

Define  the  state  transition  matrix  by: 

$(t,r)  =  T)#(t,r)  (27) 

=  I  (28) 


and  state  covariance  matrix  Ilf  =  E{s(t)2(<)^}.  As  is  well  known  [27],  the  state  covariance 
matrix  satisfies  the  following  differential  and  integral  eciuations: 

n*  =  Ftllt  +  ntFf  +  GtGf  (29) 

Et  =  $(t,0)no#(t,0)^-b  /  (.30) 

Jo 

Also,  define  the  conditional  expectation  of  the  state  at  time  ^2  given  the  states  z(ti)  and 
^(ta),  and  the  corresponding  covariance  of  the  error  as: 


Since: 


=  E{5:(/2.)k(fl),^(f3)} 

-^1!2|<i43  =  E{(2(f2)  -  3't2|<j,i3  )(^(f2)  - 

Pz{t2)\z(tl},z{t3)(^t2\^tlT  ) 


(31) 

(32) 

(33) 


20 


tlie  conditionaJ  expectation  (31)  and  error  covaidance  (32)  are  the  quantities  we  require  to 
construct  a  multiscale  representation  of  the  process  (23)  -  (26).  In  fact,  it  is  easy  to  show 
that  for  t-i  <  i'i  <  is'. 


n«j#(/2,fi)^ 

T 

nfji$(t3,ti)^ 

-1 

‘  Hh)  ■ 

^(^3,  f2)n«2 

#(f3,tl)IIij 

n^3 

Z{t3)  _ 

nq#(t2,fi)^ 

r 

IE, 

nq#(f3,fi  )^ 

-1 

'  lI/,#(t2,E)^  ' 

.  $(f3,i2)nt, 

^{tsAi  )nf, 

II«3 

^(^3:  f2)lE2 

(35) 


Equations  (34)  and  (35)  directly  provide  us  with  the  parameters  of  the  matrices  A(5), 
B(s)  and  Pq  in  the  multiscale  model  (1)  -  (3)  corresponding  to  our  mid-point  deflection 
construction.  In  particular,  let  us  identify  the  abstract  index  -s  with  a  pair  of  numbers  ( rn^ip) 
which  denote  the  scale  and  horizontal  shift  of  the  node  -s,  respectively.  The  horizontal  shift 
running  from  0  to  2'”"^  -  1,  indexes  the  nodes  at  scale  m.  For  instance,  the  root  node  is 
associated  with  the  pair  (1,0),  and  the  left  and  right  nodes  at  the  first  level  are  associated 
with  (2,0)  and  (2,1),  respectively.  With  this  notation,  the  state  at  node  .s  on  the  tree 
contains  the  values  of  the  process  r(t)  at  the  particular  three  points: 


x{s)  =  x((m,<fi)) 


z{2^t2^) 
z((2^+l}/2^) 
zi{2<p  +  2)l2^) 


(36) 


From  the  description  of  the  general  construction,  the  form  of  the  matrix  /l(.s)  in  ( 1 )  is  clear: 


/l(s)  =  A{(m,(p)) 


/  0  0 

A'l  /i2  0 

0/0 

0/0 
0  A 1  A  2 
0  0/ 


if  yp  is  even 


if  (f  is  odd 


(37) 


In  particular,  if  (p  is  even,  then  the  first  and  third  components  of  the  state  .r(s)  in  (36) 
correspond  to  the  first  and  second  components  of  a;(s7).  Thus,  the  identity  matrices  in 
(37)  for  Ip  even  simply  map  the  first  and  second  components  of  x{s^)  to  the  first  and  third 
components  of  .r(s).  In  addition,  the  mid-point  prediction  of  z({2p+  l)/2”^)  is  just  a  linear 
function  of  tlie  first  two  components  of  the  parent  .rf.sy).  wdiich  is  expressed  via  tlie  matrices 
A'l  and  A'2  in  the  second  row  of  (37).  The  matrix  /1(5)  for  p  odd  is  similar,  and  in  fact  is 
just  a  “shifted”  version  of  A{s)  for  p  even  (reflecting  the  fact  that  the  interpolation  down 
to  the  state  on  the  right  depends  on  the  last  two  components  of  x(sj)). 

The  gain  matrices  in  (37)  can  be  computed  directly  from  (34).  Using  standard  formulas 
for  the  inversion  of  a  block  2x2  matrix,  we  compute: 

Ai  =  $(/2i  tl  )  +  ^(^2!  ^l)IIij  $(t3»  )^(Ili3  —  #(  /3?  ^1  lllij  $(^3,  ti  )^)~^#(/3,  ti ) 

~n<2^(^3’^2)^(n<3  —  #(f3,ti)n<j#(t3,ti)^)“^#(f3,ti)  (38) 
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(39) 


+Iit3$(«3,^2)'-^(ni3  -  #(%,<!  )nj,#(^3,<i)^)-' 


where  h  =  2v?/2”%i2  =  i2<fi  +  l)/2™  and  =  (2ip  +  2)/2’". 
The  matrix  B{s)  in  (1)  has  the  following  block  structure: 


B(s)  =  Biim,ip)) 


0 

A3 

0 


(40) 


where  A'3  is  any  matrix  such  that  K3KJ  =  Pt2\ti,t3  and,  again,  ti  =  2vp/2“,  <2  =  (2v?+T)/2’” 
and  ta  =  (2</?+2)/2™.  The  matrix  B{s)  in  (40)  reflects  the  fact  that  ?7,o  noise  is  added  to  the 
lirst  and  third  components  of  the  state  x{s),  (which  are  simply  copied  from  the  preceding 
level),  while  noise  corresponding  to  the  error  (35)  is  added  to  the  second.  In  particidar, 
from  (3)  we  see  that  in  this  case  the  covariance  of  the  additive  noise  term  B{s)w{s}  in  (1) 
is  given  by: 


E{jB(s)i(;(s)ru^(.s)j5^(5)}  =  B{.s)B^{3) 

‘0  0  0 

=  0  -P(2v+1)/2”»|2v;/2’»,(2.v?+2)/2”'  0 

0  0  0 


(41) 


where  s  =  (in,  ip). 

Finally,  the  initial  covariance  matrix  Pq  in  (2)  is  given  by: 


Po  =  Ei 


< 

-io) 

^(0) 

1 

i 

^(0.5) 

2(0.5) 

}  (42) 

h. 

^(T) 

^(1) 

1 

llo 

$(l/2,0)no 

^(l,0)llo 


no#(i/2,0)2^  no$(i,0)2’ 

ni/2  ni/2$(  1,1/2)^ 

#(1,1/2)111/2  III 


n, 


For  instance,  if  z(t)  is  the  standard  Brownian  motion,  then  Ft 
=  t.  Thus  for  this  example  (34)  and  (35)  become: 


(43) 

0,  #(t,  r)  =  1,  and 


<2  hits 


and  we  obtain: 


A{s)  =  A{(m,(p)) 


h  -  h 
h  -  h 


^{h)  + 


h  —  h 
h  -  h 


z{t2) 


ih  -  tl)(t3  -  t2) 
h  -  h 


■  1  0  0  ■ 

1/2  1/2  0 

0  1  0 

O 

o 

0  1/2  1/2 

0  0  1 

if  ip  is  even 


if  p  is  odd 


(44) 

(45) 


(46) 
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B{s)  = 


(47) 


Po 


0 

l/2(m+l)/2 

0 

0  0  0 
0  0.5  0.5 
0  0.5  1 


(48) 


Tire  formula  (34)  for  tire  conditional  expectation  which  specifies  /l(s)  as  just  de¬ 

scribed,  also  provides  us  with  the  required  formula  for  interpolating  between  dyadic  sample 
points  at  any  level  in  our  multiscale  representation  and  hence  provides  us  with  a  direct  inter¬ 
pretation  of  this  representation  as  providing  a  sequence  of  multiresolution  approximations. 
For  example,  Brownian  motion  provides  us  with  the  linear  interpolation  formula  given  in 
(44)  and  illustrated  in  Figure  1.  This  corresponds  to  a  multiresolution  linear  spline  approxi¬ 
mation  or,  as  also  illustrated  in  Figuie  1,  as  a  non- orthogonal  multiresolution  decomposition 
using  the  so-called  “hat”  function  [57], 

As  a  second  example,  consider  the  movement  of  a  particle  whose  velocity  is  given  by  a 
Brownian  motion.  This  motion  can  be  described  using  the  following  Gauss- Markov  process: 


z{t)  = 


0  1 
0  0 


z{t)  -f- 


(49) 


In  (49),  the  first  component  of  z{t)  is  the  particle  position  and  the  second  component  is  its 
velocity.  The  state  transition  matrix  #(t,r)  and  the  state  covariance  matrix  IF  are  given 
by: 


#(f,r) 

n* 


1  t  —  T 
0  1 

i 


(50) 

(51) 


One  can  show  by  direct  computation  that  the  terms  ntj#(f2,F)^  and  #(t3,t2)n/2  in  the 
leftmost  block  matrix  on  the  right  side  of  (34)  contain  only  cubic  powers  of  t2-  Note  also 
that  the  block  matrix  in  the  middle  of  the  right  side  does  not  depend  on  #2-  Thus,  the 
interpolation  of  *(*2)  between  F  and  is  a  cubic  polynomial  in  1-2: 


Cl  -f-  C2t2  +  Cs^l  + 

C2  -f  2C3/2  +  304^2 


(52) 


where  from  (49),  the  second  component  of  Fzki.ts  Ibe  first  derivative  of  the  first. 

One  can  use  (34)  directly  to  calculate  the  values  of  the  coefficients  in  (52).  Alternatively, 
it  is  clear  from  the  definition  of  Fjpi.is  (31)  =  ^(^1)  =  ”(G)- 

These  two  constraints  provide  four  linear  equations  in  the  four  unknown  coefficients  in  (52), 
and  thus  uniquely  determine  the  interpolating  function  (52).  Note  that  the  interpolating 
polynomial  for  the  first  component  of  the  state  (the  position  of  the  particle)  has  a  continuous 
derivative  at  the  knot  locations  t  =  k/2’^,k  =  0,1, •••,2'”.  The  interpolation  of  the  first 
component  of  the  state  is  shown  in  Figures  9a  and  9b  for  the  first  two  levels  of  a  sample 
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Figure  9:  The  first  two  scales  in  a  multiscale  representation  of  a  process  which  is  equal  to 
the  second  integral  of  white  noise  are  shown.  The  representation  consists  of  samples  of  the 
process  at  dyadic  points  along  with  a  piecewise-cubic  interpolation.  Compare  these  curves 
with  the  lower  two  graphs  of  Figure  1,  which  depict  the  piecewise  linear  interpolation  of 
the  first  integral  of  white  noise. 
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path  of  the  multiscaJe  realization.  Figure  9a  is  the  cubic  interpolation  of  <r(0),^(0.5)  and 
2(1),  while  Figure  9b  is  the  cubic  interpolation  of  ^(0),  0.25),  ^:{0.5),  *(0.75)  and  ^(l). 

Finally,  while  we  have  focused  on  the  construction  of  multiscale  models  for  continuous- 
time  Gauss-Markov  processes,  an  exactly  analogous  set  of  calculations  can  be  performed 
for  the  discrete-time  process: 


=  Ftzit)  (.53) 

Also,  as  discussed  in  Section  3.2,  in  this  case  we  can  either  construct  models  with  three- 
point  state  vectors  (e.g.  [*:(0),  z{Tl2),  *(T)]  or  four-point  state  vectors  [*(1),  z{T/2),  z{T/2+ 
1 ),  z(T )] ).  The  former  of  these  is  exactly  analogous  to  what  we  have  done  here  in  cojitinuous- 
time  and  has,  as  we  have  indicated,  a  high  degree  of  redundancy,  wdiile  the  latter  does  not. 
We  defer  explicit  discussion  and  illustration  of  such  non-redundant  i-epresentations  until 
Section  4  where  we  describe  them  in  the  context  of  modeling  2-D  MRF’s. 


3.3.2  Discrete-State  Processes 


Next  consider  a  general  finite-state  Markov  process  2( t)  G  {1 ,2,  ■  •  ■  L]  defined  over  a  discrete 
interval  /  G  {0, 1,  ■  •  • ,  2’}.  The  probability  structure  of  the  process  is  completely  determined 
by  the  initial  conditions: 

Pr[*(0)  =  k]  (54) 

for  k  G  {1,2,  “  ■£}  and  by  the  one-step  transition  probabilities: 


Fij  =  Pr[z(t)  =  ilz(t-l)  =  J] 
We  define  the  one-step  transition  matrix: 

•  •  •  Fi,l  ■ 

•  •  •  F2,L 

•  •  •  ^L,L  . 


A,1  Fi,2 

F2,1  P2,2 

Pl,i  Pl,2 


(55) 


(56) 


Note  that  the  multistep  transition  probabilities  are  given  by  powers  of  the  matrix^'^  P: 


Pr[^(t  4-  r)  =  i\z{t)  =  j]  =  [P^]i,j 


(57) 


Using  (57)  and  Bayes’  rule  it  is  straightforward  to  calculate  that  for  <  t2  <  #3: 

PrN^2)  =  Mh)  =  u  i 

These  conditional  probabilities,  in  addition  to  the  probability  function  required  for  the  state 
at  the  root  node  of  the  tree,  namely 

FiiziQ)  =  £z(T/2)  =  j,z{T)  =  k]  =  [F^/']fcj[2’^/"k,:Prb(0)  = /]  (-59) 

stands  for  the  {i,j)  element  of  the  matrix  A. 
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Figure  10:  A  sample  patli  of  a  discrete-state  Markov  process  is  shown.  A  multiscale  repre¬ 
sentation  of  this  8  state  process  is  developed  in  Section  3.3.2. 


allow  us  to  construct  the  multiscale  representation  of  the  process.  Note  that  (58)  is  the 
counterpart  of  the  conditional  probability  equations  for  Gauss-Markov  processes  given  in 
(33)  -  (35),  and  that  the  pdf  for  the  state  at  the  root  node  (59)  is  the  counterpart  of  the 
initial  covariance  matrix  (43). 

One  special  case  of  this  process  is  the  following: 


(60) 

p. .  _ 

(61) 

Piiz{0)-ij  =  l/£,  i  = 

(62) 

A  sample  path  of  such  a  process  is  shown  in  Figure  10  for  L  =  8  and  /t  =  0.97.  Neighboring 
states  of  this  process  tend  to  be  the  same,  and  when  the  process  does  change  state,  no 
particular  change  is  preferred.  Thus,  this  model  would  seem  to  be  a  natural  one  to  use  in 
segmentation  applications  and  can  in  fact  be  viewed  as  an  alternative  to  the  1-D  multiscale 
model  (6)  introduced  in  [7,  8,  9].  As  noted  in  Section  2.2,  the  model  in  (6)  does  tiot  in 
general  produce  a  Markov  chain  or  reciprocal  process  at  the  finest  level.  On  the  other  hand 
(60)  -  (62)  is  a  Markov  model  and  for  this  process: 


(1-f  (T-l)#)/X  if?:=:i 
(  l-d^  )/£  if  i  ^  j 


(63) 
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where: 


i)  =  {Ln-l)/iL-l) 


(64) 


The  conditional  probability  function  (63)  can  be  verified  by  noting  that,  given  (60),  (61), 
the  one-step  transition  matrix  (56) is  circtilant,  and  thus  diagonalized  by  the  LxL  Discrete 
Fourier  Transform  matrix.  Alternatively,  by  writing  the  one-step  transition  matrix  for  this 
example  as  the  sum  of  a  multiple  of  the  identity,  plus  a  matrix  with  constant  entries,  one  can 
directly  calculate  the  eigenvalues  and  eigenvectors,  which  again  can  be  used  to  diagonalize 
(56).  ‘ 

Using  (58),  for  this  example  we  can  write  down  the  transition  probabilities  for  the  mid¬ 
point  deflection  model.  In  particular,  assuming  that  T  is  a  power  of  two,  we  can  associate 
the  state  at  node  s  with  the  following  values  of  the  process: 


x(s)  =  .x-((m,  <fi)) 


z{2<pTfT^) 
z({2<pT  +  T)/2"^) 
z({2^T  +  2T)l2^'^) 


(65) 


where,  as  in  (36),  the  pair  of  numbers  (m,  </:)  denote  the  scale  and  horizontal  shift  of  the 
node  s,  respectively.  Thus,  to  generate  the  state  at  node  s,  given  the  state  at  the  parent 
node  .S7,  we  require  the  following  conditional  pdf: 


C1C1/C2 

if  i 

=Z 

j  =  k 

Pr[^ 

2^T  +  T  2^T 

'  Orn  ■  1  V  fyxfi  ' 

..  2<^T-f2r 

)  =  k]=< 

C1C1/C2 

C1C1/C2 

if  i 
if  i 

j  =  k 
j  #  k 

(66) 

C1C1/C2 

if  i 

= 

k  7^  j 

C1C1/C2 

if  i 

1  J' 

k  distinct 

where; 

Cl 

=  (1  +  (T- 

l)^r/2- 

')/L 

(67) 

C2 

=  (l  +  (i- 

l)^r/2- 

^)/L 

(68) 

Cl 

=  (1  - 

(69) 

C2 

(70) 

To  gain  additional  insight  concerning  the  structure  of  our  multiscale  models,  consider 
the  particular  example  of  a  stationary  two-state  binary  process  with  one-step  transition 
inatrix  and  initial  state  probabilities  equal  to: 


P  = 


Pr[^(0)  =  i]  = 


1  -  //,  1] 

M  1  -  _ 

rj/iV  +  ft)  if  (  =  1 

+  ifj  =  2 


For  this  pi’ocess  one  can  show  that: 


pk 


1 


r]  +  n  f  rj  +  ?/( 1  -  ?/  - 

+  IJ,{1  -  t]  -  fj,)^  H  +  T](l  -  7]  - 


(71) 

(72) 


(73) 
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and  tliiis  using  (73)  with  (58),  (59)  one  can  build  multiscale  representations  for  the  class 
of  stationary  binary  Markov  processes.  While  we  have  focused  in  this  example  and  in  the 
previous  one  on  processes  with  stationary  state  transition  probabilities,  it  is  straightforward 
to  apply  this  construction  to  the  representation  of  non-stationary  discrete-state  Markov 
processes  as  well,  simply  by  choosing  the  conditional  probability  functions  required  in  the 
multiscale  representation  correctly. 

Moreover,  the  mid-point  deflection  structure  can  also  be  used  to  generate  non-Markov 
processes  on  the  tree.  For  instance,  consider  the  following  binary  mid-point  “selection” 
process  deftned  over  t  E  {0, 1,  •  •  -2^}  [59]: 

Pr[r(0)  =  ?:,5(2^-\)  =  i,^(2^)  =  A:]  =  1/8  foralli,j,A:  (74) 

(  1.1  if  i  =  j  =  k 

Pr[2(i2)  =  i\zitx)  =  i,  z{h)  =  k]  =  <  1  -  /tt  if  i  i  and  j  =  k  (75) 

(  0.5  \i  j  ^  k 

where  i,j,k  E  {1,2}  and  where  ti,t2i^3  are  any  3-tuple  of  dyadic  points  corresponding 
to  one  of  the  state  vectors  in  the  multiscale  representation.  At  the  coarsest  “scale”  of 
this  process,  the  three  components  of  the  state  vector  xq  are  independent  and  identically 
distributed  random  variables,  each  equally  likely  to  be  1  or  2.  It  is  ea,sy  to  show  that 
the  process  resulting  from  this  construction  is  not  Markov  in  general,  and  thus  we  can 
coixclude  that  the  set  of  binary  stochastic  processes  which  can  be  constructed  within  the 
mid-point  deflection  framework  is  strictly  larger  than  the  class  of  binary  Markov  processes 
over  intervals. 

In  fact,  a  bit  of  thought  shows  that  the  class  of  processes  realizable  by  mtdtiscale  models 
is  quite  a  bit  larger  than  the  class  of  Markov  chains.  Indeed,  any  binary  stochastic  process 
defined  over  t  E  {0, 1,  •  •  ■  ,2^}  when  represented  via  mid-point  deflection  has  a  probability 
structure  which  is  determined  by  4(2^  -  1)  parameters,  corresponding  to  the  required 
conditional  probability  functions.  In  particular,  the  conditional  probabilities  Pr[r(t2)  = 
=  j,z(t3)  =  k]  for  specific  choices  of  ti  <  t.2  <  h  are  uniquely  determined  by  the 
four  parameters: 


Pr[a(t2)  =  l|a(ti)  =  l,^(/3)  =  1]  =  Aj  (76) 

Fr[zit2)  =  l\z{ti)  =  l,z(t3)  =  2]  =  X2  (77) 

PT[zit2)  =  l\z(ti)  =  2,z{t3)  =  l]  =  As  (78) 

Pr[a(/.2)  =  l|a(ti)  =  2,^(t3)  =  2]  =  A4  (79) 


Since  the  process  is  represented  with  an  iV  level  mnltisrale  process.  Ihere  are  2^  —  2  of  these 
conditional  densities  which  must  be  specified,  corresponding  to  each  of  the  nodes  except  the 
root  node.  However,  the  probability  function  for  the  state  at  the  root  node  also  requires  four 
parameters,  and  thus  the  total  number  of  parameters  to  be  specified  is  4(2^^— 1 ).  In  contrast, 
a  non-stationary  binary  Markov  process  defined  over  the  time  interval  t  E  {0, 1,  •  •  • ,  2^^} 
requires  at  most  1  -f  2  X  2^  parameters  (one  corresponding  to  the  initial  probability,  and 
2  for  each  transition  from  t  to  t  -|-  1,  for  t  =  0, 1,  •  •  •2’^  -  1).  Since  each  of  the  parameters 
in  each  of  these  models  is  a  probability,  i.e.  a  number  in  the  interval  [0,1],  we  see  that  the 
set  of  processes  arising  from  N -level  multiscale  models  is  in  one-to-one  correspondence  with 
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the  4(2^-  l)-dimensional  unit  cube,  while  the  set  of  non-stationary  Markov  chains  over  the 
same  length  interval  (2^  +  1)  corresponds  to  the  2(2^  +  l)-dimensional  unit  cube.  Thus, 
for  N  >  I,  Markov  processes  constitute  only  a  “thin”  subset  of  the  entire  class  of  binary 
tree  processes. 

4  Represention  of  2-D  Markov  Random  Fields 

In  this  section  we  first  review  a  few  of  the  properties  of  Markov  random  fields  and  then 
describe  how  Markov  random  fields  can  be  represented  exactly  with  our  multiscale  modeling 
structure.  Next,  we  introduce  a  family  of  approximate  representations  for  Gaussian  MRF’s 
employing  1-D  wavelet  transforms. 

4.1  2-D  Markov  Random  Fields 

Markov  random  fields  (MRF’s)  are  a  multidimensional  generalization  of  1-D  reciprocal 
processes.  A  continuous  space  stochastic  process  £  7^"  is  said  to  be  a  Markov 

random  field  if  the  probability  distribution  of  the  process  in  the  interior  of  any  closed  set 
D  is  independent  of  the  process  outside,  conditioned  on  the  values  of  z(t)  on  the  boundary 
r  of  That  is,  for  t  6  D  \  F: 

Pz(t)ls(r),re(^i\r)4^tl^T,T  6  (fl  \  rf)  = 

^  T)  (80) 

wlxere  the  notation  D  \  F  denotes  the  set  of  elements  in  D  which  are  not  in  F  (in  this  case, 
the  interior  of  Q).  The  definition  for  Markov  random  fields  on  discrete  lattices  requires 
the  specification  of  the  notion  of  the  “boundary”  of  a  set  in  2"  [64,  29].  Typically,  this 
is  accomplished  through  the  specification  of  a  neighborhood  system.  The  idea  is  that  the 
probability  distribution  of  z(t),  conditioned  on  a  set  {*(r),r  G  Dt]  in  the  neighborhood, 
Dt,  of  t,  is  independent  of  the  process  outside  the  neighborhood: 

Pz{t)\z{T),T&Z’^\{tyi  ^  G  2  \  "{0}  ~ 

Pz{t)\z{T},TEDti^f\^Tl'^  G  (81) 

In  this  paper,  we  focus  on  2-D  MRF’s,  i.e.  where  t  G  2'^,  and  in  this  context  there  is 
a  hierarchical  sequence  of  neighborhoods  frequently  used  in  image  processing  applications 
[13].  The  fir.st  order  neighborhood  of  a  lattice  point  consists  of  its  four  nearest  neighbors 
(in  the  Manhattan  metric),  and  the  second-order  neighborhood  consists  of  its  eight  nearest 
neighbors.  The  sequence  of  neighborhoods  up  to  order  seven  is  illustrated  in  Figure  11. 

A  given  neighborhood  system  implicitly  determines  the  boundary  set  of  any  particular 
region.  In  particular,  given  the  neighborhood  system  G  Z'^.  the  boundary  F  of  a  subset 
0  of  2^  is  given  by  the  set  of  points  which  are  neighbors  of  elements  in  fl,  but  not  elements 
of  Q: 


F  =  {r|r  G  A,/ G  0}  \  fi  (82) 

The  conditional  distribution  and  neighborhood  structure  cannot  be  chosen  arbitrarily 
if  one  is  to  obtain  a  consistent  joint  distribution  for  the  elements  of  the  random  field. 
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Figure  11:  The  first  through  seventh-order  neighborhoods  of  lattice  site  t  are  shown.  The 
first-order  neighborhood  consists  of  just  the  two  vertical  and  two  horizontal  nearest  neigh¬ 
bors. 


First,  the  neighborhood  system  must  have  the  properties  that  t  ^  Dt  and  (2)  if  /  G  Hr 
then  T  G  Df  Second,  the  conditional  distribution  functions  must  satisfy  the  consistency 
conditions  given  by  the  Hammersley-Clilford  theorem  [4].  For  detailed  accounts  of  these 
issues  and  MRF’s  in  general,  we  refer  the  reader  to  a  few  of  the  widely  referenced  papers 
in  the  field  [30,  54,  56,  4,  64,  36,  33,  29]. 

4.2  Exact  Multiscale  Representations  of  2-D  Markov  Random  Fields 

The  representations  of  1-D  reciprocal  and  Markov  processes  in  Section  3  relied  on  the 
conditional  independence  of  regions  inside  and  outside  a  boundary  set,  and  we  use  the  same 
idea  here  to  represent  Markov  random  fields  on  a  stiuare  lattice.  The  multiscale  model  is 
identical  to  that  used  in  the  1-D  case,  except  that  it  is  defined  on  a  quadtree  instead  of 
a  dyadic  tree.  That  is,  we  consider  multiscale  models  exactly  as  in  (1)  -  (3)  but  where  $ 
denotes  a  node  on  the  ciuadtree  depicted  in  Figure  12  and  7  is  a  four-to-one  operator,  i.e. 
each  point  is  the  parent  of  four  descendant  points  at  the  next  level. 

Quadtree  structures  arise  naturally  in  multiresolution  image  processing  contexts.  For 
instance,  successive  filtering  and  decimation  operations  lead  to  images  defined  on  such  a 
hierarchy  of  grids  in  the  Laplacian  pyramid  coding  algorithm  of  Burt  and  Adelson  [11]  and 
in  the  closely  related  wavelet  transform  decomposition  of  images  [44].  .41so.  the  multigrid 
approaches  to  low  level  vision  problems  discussed  by  Terzopoulos  [60]  involve  relaxation  on 
a  similar  seciuence  of  grids.  In  addition,  Wilson  and  his  colleagues  [21]  have  introduced 
and  used  a  particular  simple  example  of  the  model  in  (1)  -  (3)  with  the  state  defined  on  a 
quadtree  for  image  processing. 

In  fact,  it  is  interesting  to  note  that  both  in  [21]  and  in  the  segmentation  studies  in 
[7,  8,  9]  a  similar  deficiency  in  the  quadtree  models  was  discussed.  In  particular,  processing 
based  on  these  models  led  to  blocky  results  which  could  differ  significantly  for  different 
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Figure  12:  The  quadtree  structure  shown  is  used  for  the  multiscale  representations  of 
Markov  random  fields  (MRF’s).  Each  node  of  the  quadtree  has  four  offspring,  denoted 
saNWf'^f^NE^sasEi^f^sw-  Again,  the  parent  of  node  s  is  denoted  37,  and  in  this  case  7  is 
a  four-to-one  shift  operator. 


positionings  of  the  tree  with  respect  to  the  finest  scale  Z'^  lattice,  and  this  fact  led  to  the 
need  for  more  complex  processing  structures.  In  particular,  in  [21]  it  was  necessary  to 
perform  processing  several  times  with  different  tree  positions  and  to  average  the  results, 
while  in  [7,  8,  9]  the  tree  was  replaced  by  a  more  connected  lattice  (so  that  any  two  nodes 
have  both  common  ancestors  and  common  descendants).  The  latter  modification,  however, 
destroys  the  partially-ordered  Markovian  structure  which  the  dyadic  and  quadtree  processes 
possess  and  which  leads  to  highly  parallelizable  and  scale- ?’ecM.r.sfue,  rather  than  iterative, 
algorithms.  In  this  section,  we  show  that  one  can  completely  avoid  the  apparent  problems 
in  using  quadtree  models  by  demonstrating  that  one  can  model  any  MRF  exactly  using  such 
a  model. 

Consider  a  2-D  MRF  z{t)  defined  on  a.  (2^  -|- 1)  x  (2^  -h  1)  lattice.  The  construction  of 
reciprocal  processes  in  one-dimension  started  with  the  values  of  the  process  at  the  initial, 
middle  and  end  points  of  an  interval.  In  two  dimensions,  the  analogous  top  level  description 
consists  of  the  values  of  the  MRF  around  the  outer  boundary  of  the  lattice  and  along  the 
vertical  and  horizontal  “mid-lines”  which  divide  the  lattice  into  four  quadrants  of  equal  size. 
For  instance,  on  a  17x  17  lattice,  the  state  vector  Xq  at  the  root  node  of  the  quadtree  contains 
the  values  of  the  MRF  at  the  shaded  boundary  and  mid-line  points  shown  in  Figure  13. 
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o  Boundary  points 


o  Mid-line  points 


Figure  13:  The  state  vector  at  the  root  node  in  the  MRF  multiscale  representation  consists 
of  the  MRF  values  at  the  boundary  and  “mid-line”  points,  shown  in  the  shaded  region  here 
for  a  17  X  17  lattice.  To  construct  a  sample  path  of  the  MRF  using  the  “mid-line”  deflection 
construction,  we  start  by  choosing  a  sample  from  the  joint  distribution  of  the  values  in  the 
root  node  state. 

The  boundary  points  are  denoted  with  O  and  o  symbols,  respectively.  In  general,  the  state 
at  the  root  node  is  a  (6  X  2^  -  3)-dimensional  vector  (given  some  ordering  of  the  boundary 
and  mid-line  lattice  points).  To  construct  a  sample  path  of  the  MRF,  we  Imgin  by  choosing 
a  sample  from  the  joint  pdf  of  the  MRF  values  defined  on  the  boundary  and  mid-line  set. 
This  is  the  2-D  counterpart  to  choosing  a  sample  from  the  pdf  in  (20)  when  constructing  a 
1-D  reciprocal  process. 

In  the  1-D  case,  transitions  from  the  first  to  second  level  consisted  of  obtaining  a  sample 
from  the  conditional  distribution  of  the  state  at  the  mid-points  of  the  left  and  right  half¬ 
intervals.  In  two  dimensions,  we  predict  the  set  of  values  at  the  mid-lines  in  each  of  the 


four  quadrants.  The  components  of  the  four  state  vectors  at  the  second  level  are  illustrated 
in  Figure  14  for  the  17  X  17  MRF.  The  points  corresponding  to  the  state  in  the  north-west 
corner  are  shaded,  and  correspond  to  a  scaled  and  shifted  version  of  the  points  at  the  top 
level.  The  boundary  points  of  the  north-west  state  are  denoted  with  open  and  blackened 
diamond  symbols  and  the  new  mid-line  points  are  denoted  with  open  circles.  Note  that  the 
four  states  at  the  second  level  share  the  black  diamond  mid-line  poiirts  of  the  state  at  the 
first  level.  This  is  analogous  to  the  1-D  construction  in  which  the  mid-point  at  the  first  level 
corresponds  to  an  end  point  in  both  states  at  the  second  level  (cf.  Figure  3).  In  general,  the 
state  vectors  at  the  four  nodes  at  the  second  level  each  specify  the  MRF  at  6  X  2^  —  3 

lattice  points.  Each  of  these  states  at  the  second  level  consists  of  points  carried  down  from 
the  root  node  (namely  the  diamond  boundaries  of  each  of  the  quadrants  in  Figure  14 )  as 
well  as  new  mid-line  points  within  each  quadrant  (the  open  circles  in  Figure  14).  These 
mid-line  values  are  chosen  as  samples  from  their  joint  conditional  distribution,  given  the 
state  at  the  root  node.  The  key  point  here  is  that  given  the  values  of  the  field  around  the 
boundary  of  each  quadrant,  the  values  of  the  field  along  the  mid-lines  of  that  quadrant 
are  independent  of  the  values  outside  this  ciuadrant.  Said  another  way,  the  four  states  at 
the  second  level  of  the  tree  are  conditionally  indejiendent  given  the  values  of  the  MRF  on 
their  respective  boundaries,  i.e.  given  precisely  that  information  captured  in  the  state  at 
the  first  level.  Thus,  the  values  along  the  new  mid-lines  at  the  second  level  can  be  chosen 
independently  and  in  parallel,  in  analogy  to  the  way  the  two  mid-points  in  (21),  (22)  are 
chosen. 

Now,  we  can  iterate  the  construction  by  defining  the  states  at  successive  levels  to  be 
the  values  of  the  MRF  at  boundary  and  mid-line  points  of  successively  smaller  subregions. 
Indeed,  by  subdividing  each  quadrant  in  the  same  way  as  we  did  in  going  from  the  fust  level 
to  the  second,  at  the  level  the  state  vectors  each  contain  the  values  of  the  MRF  at 
6  X  —  3  boundary  and  mid-line  points.  Note  that  the  dimension  of  the  state  varies 

from  level  to  level,  reflecting  the  obvious  fact  that  the  number  of  points  in  the  boundary  of 
a  2-D  region  depends  on  the  size  of  the  region.  The  multiscale  representation  has  N  levels, 
and  each  of  the  4^“^  states  at  level  N  represent  9  values  in  a  3  x  3  square.  Because  of  the 
Markov  property,  at  each  level  the  states  are  conditionally  independent,  given  their  parent 
state  at  the  next  higher  level.  Thus,  the  MRF  can  be  thought  of  precisely  as  a  multiscale 
stochastic  process,  and,  in  the  Gaussian  case,  this  leads  to  models  exactly  as  in  (1)  --  (3). 

As  in  the  1-D  case,  there  are  several  comments  to  make.  First,  we  have  described  a 
construction  in  which  the  lattice  is  square.  If  the  MRF  is  defined  over  a  non-square  lattice, 
then  the  same  basic  approach  will  work.  In  particular,  all  we  require  is  some  sequence  of 
subregions  whose  boundaries  eventually  become  dense  in  the  set  of  lattice  points.  Second, 
just  as  our  1-D  multiscale  model  has  a  natural  interpretation  in  terms  of  decimation  — 
e.g.  if  the  points  on  the  finest  scale  correspond  to  integers,  i.e.  to  Z.  then  at  the  next  most 
fine  scale  they  correspond  to  even  integers,  i.e.  22  —  so  does  our  2-D  model,  although  it 
differs  from  the  usual  notion  of  decimation  in  2-D.  Specifically,  if  the  points  on  the  finest 
scale  correspond  to  2^  =  2  x  2,  then  the  usual  notion  of  decimation  would  be  22  x  22. 
In  contrast,  the  notion  of  decimation  associated  with  our  multiscale  models  yields  the  set 
(22  X  2)(J(2  X  22)  at  the  next  finest  scale. 

Third,  note  that  the  representation  we  have  defined  is  redundant.  In  the  1-D  case  this 
redundancy  was  removed  by  predicting  two  mid-points  at  each  level  instead  of  one  (cf. 
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o  Boundary  points  at  both  first  and  second  levels 

♦  Boundary  points  at  the  second  level  and  mid-line 
points  at  the  first  level 

o  New  (second  level)  mid-line  points 


Figure  14:  The  components  of  the  four  state  vectors  at  the  second  level  of  the  tree  are  scaled 
and  shifted  versions  of  the  components  of  the  state  at  the  root  node.  For  instance,  the  state 
corresponding  to  the  north-west  corner  at  the  second  level  of  a  representation  for  an  MRF 
defined  on  a  17  x  17  lattice  consists  of  the  values  of  the  process  at  the  shaded  points.  The 
values  of  the  MRF  at  the  boundary  points  in  these  second  level  states  are  mapped  dowm 
from  the  root  node  state,  and  the  values  at  the  new  mid-lines  in  each  of  the  four  quadrants 
are  chosen  independently.  In  particular,  the  new  mid-fine  values  in  any  given  quadrant  are 
independent  of  values  of  the  MRF  outside  that  quadrant,  given  the  boundary.  Thus,  in 
the  construction  of  a  sample  path,  we  can  choose  values  along  each  of  the  four  sets  of  new 
mid-lines  independently  and  in  parallel.  This  process  can  then  be  iterated,  by  defining  the 
states  of  the  multiscale  process  at  lower  levels  in  the  quadtree  with  respect  to  successively 
smaller  subdomains,  and  constructing  the  process  (along  boundary  and  mid-line  points) 
independently  within  each  subdomain. 
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Figure  6).  Likewise,  in  two  dimensions  we  can  eliminate  the  redundancy  by  predicting  ttoo 
mid-lines  at  each  level  instead  of  one.  For  instance,  the  state  vector  at  the  top  level  of  the 
quadtree  for  an  MRF  defined  on  a  16  x  16  lattice  would  consist  of  the  values  of  the  MRF 
at  the  shaded  points  shown  in  Figure  15,  and  the  four  state  vectors  at  the  next  level  would 
consist  of  the  values  of  the  process  at  the  points  shaded  in  Figure  16.  Iir  the  2-D  case  this 
latter  representation  is  important  in  the  context  of  defining  approximate  models,  as  will  be 
seen  in  the  next  section,  and  illustrated  in  Section  5. 

Fourth,  higher-order  models  can  be  easily  accommodated.  For  example,  if  the  Markov 
random  field  has  up  to  a  fifth-order  neighborhood  (cf.  Figure  11),  the  field  can  be  repre¬ 
sented  by  taking  as  state  the  values  of  the  process  along  boundaries  and  mid-lines  of  “width” 
two.  In  general,  Jieighborhoods  of  any  order  can  be  handled  by  increasing  the  boundary 
and  mid-line  widths  appropriately. 

Fifth,  in  the  case  of  Gaussian  MRF’s,  the  prediction  from  one  level  to  the  next  in  our 
representation,  as  captured  by  the  term  ).r(s7)  in  ( 1 ),  is  simply  the  result  of  a  2-D  version 
of  the  interpolation  formula  displayed  in  (31)  for  general  1-D  Gauss-Markov  processes  and 
in  particular  in  (44)  for  Brownian  motion.  As  in  1-D,  this  interpolation  formula  can  also 
be  thought  of  as  providing  a  multiresolution  approximation  to  a  process  using  interpolation 
functions  specifically  tailored  to  the  process  under  consideration.  To  see  this  more  clearly, 
note  that  the  linear  spline  interpolation  formula  for  Brownian  motion  given  values  at  two 
points  ~(0)  =  Zo  and  z{T)  =  Zt  is  simply  the  solution  to  the  second-order  differential 
equation : 


^^o,T  - 


0 


(S3) 


Similarly  the  interpolation  of  the  first  component  of  the  second-order  process  (49)  is  given 
by  the  solution  of: 


0 


(84) 


given  4:(0),  i(0),  ^(r)  and  z{T}.  The  2-D  example  analogous  to  the  linear  spline  model  for 
Brownian  motion  is  Laplace’s  equation: 


=  0 


(85) 


given  values  of  5  on  the  boundary  of  a  square  region,  while  the  counterpart  to  (84),  cor¬ 
responding  to  a  second-order  model,  would  be  the  solution  of  a  homogeneous  biharmonic 
equation: 

=  0  (86) 


given  boundary  values  and  normal  derivatives  along  the  boundary.  More  generally,  the 
particular  interpolation  formula  for  a  1-D  or  2-D  process  is  given  by  the  solution  of  a  specific 
homogeneous  differential  (or  partial  differential)  equation  determined  by  the  covariance 
structure  of  the  process  (see,  for  example  [63,  39],  for  related  discussions). 

Finally,  we  note  that  for  domains  of  substantial  size,  the  representations  may  be  of 
prohibitively  large  dimension.  This  issue  is  addressed  for  Gaussian  MRF’s  in  the  next 
section,  where  we  introduce  a  family  of  low-dimensional  approximate  representations  based 
on  one-dimensional  wavelet  transforms  of  the  MRF  along  1-D  boundaries. 
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Figure  15:  The  state  at  the  root  node  in  a  non-redimdant  exact  multiscale  representation 
of  an  MRF  defined  on  a  16  X  16  lattice  consists  of  the  values  of  the  process  at  the  sliaded 
points.  The  redundancy  in  the  exact  representation  is  eliminated  by  generating  the  values 
of  the  process  along  two  mid-lines  instead  of  one.  The  figure  also  illustrates  the  sets 
and  the  sequences  jda,i,i{k)  defined  in  the  context  of  approximate  representations  in  Section 
4.3.  The  are  1-D  sequences  corresponding  to  values  of  the  MRF  along  boundaries 

of  square  subdomains  (which,  at  the  first  level,  are  the  white  areas  in  the  figure).  These 
sequences  overlap  at  the  corner  points  of  Ijoundaries.  In  the  figure,  this  is  represented  by 
putting  two  symbols  at  the  same  lattice  point,  e.g,  y  t>  di  the  upper  left  corner. 
The  approximate  representations  take  as  the  state  subsets  of  the  coefiftcients  in  1-D  wavelet 
expansions  of  the  /3s,i,j{k)  sequences. 
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Figure  16:  The  four  states  at  the  second  level  of  the  tree  in  a  non-redundant  exact  multiscale 
representation  are  scaled  and  shifted  versions  of  the  state  at  the  root  node,  and  are  shown 
here  for  an  MRF  defined  on  a  16  X  16  lattice.  The  state  in  the  north-west  corner  contains 
the  values  of  the  process  at  the  shaded  points  in  the  north- w'est  8x8  quadrant.  With  the 
node  .s  corresponding  to  this  north-west  corner  state,  the  sets  F^.j  and  sequences  lh,NW,j 
are  illustrated.  Note  again  that  the  sequences  j3s,i,j  overlap. 


37 


4.3  Approximate  Multiscale  Representations  of  2-D  Gaussian  Markov 
Random  Fields 

In  this  section  we  propose  a  famiiy  of  approximate  representations  for  Gaussian  MRF’s  that 
provide  iow- dimensional  aiternatives  to  the  exact  muitiscaie  representations.  The  basic  idea 
behind  the  approximate  representations  is  to  take  as  the  state  not  boundaries  of  regions,  but 
rather  some  reduced-order  representation  of  them.  Conceptualiy,  we  wouid  iike  to  retain 
oniy  those  components  of  the  boundary  that  are  rec^uired  to  maintain  neariy  compiete 
conditional  independence  of  regions.  In  general,  exact  conditional  independence  wiU  be  lost 
unless  the  entire  boundary  is  kept,  but  as  we  discuss  and  illustrate  here  and  in  the  next 
section,  in  many  cases  only  a  small  amount  of  information  needs  to  be  retained  in  order 
to  obtain  adequate  representations  of  the  important  statistical  and  qualitative  features  of 
a  Gaussian  MRF. 

The  basis  for  our  approximation  methodology  is  a  change  of  coordinates  in  representing 
the  values  of  MRF’s  along  1-D  boundaries.  A  family  of  models  can  then  be  generated  by 
making  different  choices  for  the  set  of  coordinates  to  be  retained  and  those  to  be  discarded 
at  each  level  of  the  multiscale  representation.  These  models  range  from  being  exact  (if 
all  coordinates  are  retained)  to  increasingly  approximate  and  simple  as  fewer  and  fewer 
coefficients  are  retained.  While  one  can  imagine  using  a  wide  variety  of  different  coordinate 
transformations,  including  a  1-D  Fourier  transform  or  a  Karhunen-Loeve  expansion,  we 
focus  here  on  a  choice  that  is  particularly  well  matched  to  the  self-similar,  multiresolution 
nature  of  our  exact  representation.  Specifically,  as  illustrated  in  Figures  15  and  16,  as  we 
proceed  from  one  level  to  the  next  finer  level  in  our  multiscale  representation  of  a  Gaussian 
MRF,  the  boundaries  at  the  coarser  levels  are  essentially  halved  in  length  (and  new,  shorter 
boundaries  are  added  as  well).  This  self-similar  structure  suggests  representing  the  values 
of  the  MRF  along  such  boundaries  in  terms  of  wavelet  bases.  In  this  case,  if  at  each  level 
we  keep  oidy  the  wavelet  coefficients  up  to  a  particular  level  of  resolution,  then  at  each  level 
in  our  representation  we  are  actually  adding  only  one  new  level  of  detail. 

The  approximation  that  is  made  in  keeping  only  coarse  wavelet  approximations  along 
boundaries  actually  has  two  parts.  The  first  is  that  we  assume  that  the  new  detail  to  be 
added  along  each  boundary  from  level  to  level  is  independent  of  the  previously  generated 
coarse  approximatioixs.  The  second  is  that  we  neglect  the  residual  correlation  between  MRF 
values  in  neighboring  subregions  when  we  are  given  only  a  coarse  approximation,  rather 
than  complete  knowledge  of,  the  values  along  their  common  boundary.  The  former  of  these 
approximations,  namely  the  scale-to-scale  decorrelation  capability  of  wavelet  transforms, 
has  already  been  studied  and  well  documented  in  several  papers  on  1-D  stochastic  processes 
[19,  61].  The  latter  of  these,  which  deals  explicitly  with  the  full  statistical  structure  of  2- 
D  MRF’s,  has  not,  to  our  knowledge,  been  investigated  previously  (in  fact,  the  use  of 
1-D  wavelets  for  2-D  random  fields  is,  we  believe,  a  completely  new  idea).  As  the  results 
presented  here  illustrate,  this  appears  to  be  an  extremely  effective  method  for  many  MRF's. 

The  approximate  models  are  based  on  the  non-redundant  exact  representations  for 
MRF’s  described  in  the  previous  section.  The  states  at  the  first  and  second  levels  of  this 
representation  for  an  MRF  defined  on  a  16  x  16  lattice  are  shown  in  Figures  1-5  and  16. 
The  root  node  state  in  Figure  15  contains  the  values  of  the  MRF  at  112  points.  More 
generally,  in  a  multiscale  representation  of  an  MRF  defined  on  a  2^  x  2^  lattice,  a  state  at 
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the  level  in  this  exact  representation  represents  the  values  of  the  MR.F  at  —  1) 

points.  We  denote  this  set  of  points  as  and  we  view  it  as  the  union  of  four  mutually 
exclusive  subsets.  In  particular,  consider  again  the  112  points  associated  with  the  root  node 
state  in  Figure  15.  We  can  view  these  as  four  sets  of  28  points,  each  of  which  corresponds 
to  the  boundary  of  one  8x8  quadrant.  In  general,  we  can  divide  the  set  F^  into  four 
sets  of  4(2^”"d»)  —  1)  points  in  a  similar  fashion,  and  we  denote  these  mutually  exclusive 
subsets  as  Fs,i, i  6  {NW,N E^SE^SW}^  where  the  subscripts  refer  to  the  spatial  location 
of  the  subset.  For  instance,  with  s  =  0  corresponding  to  the  root  node,  the  four  subsets 
Fo,i,*  6  (NW,N E,SE,SW}  are  illustrated  in  Figure  15  with  the  symbols: 


Fo,ivvv  V? 

Fo.ate  0 

Fo,5je;  ^  0 

Fo,5T'F  ^  ° 

<1,  t>.  A,  and  combinations  of  these. 

(87) 

(88) 

(89) 

(90) 

(91) 

Next,  we  interpret  the  set  of  values  {2(f),  t  G  F^.i}  for  each  of  these  quadrant  boundaries, 
as  fonr  1-D  sequences  of  length  corresponding  to  each  of  the  sides  of  the  quadrant 

boundary.  Thus,  there  are  a  total  of  sixteen  1-D  boundary  sequences  associated  with  the 
set  Fj,,  and  we  denote  these  as:  €  {NW,NE,SE,SW},j  G  {hu^  hi,  vl,vr},  where 

the  latter  four  subscripts  refer  to  the  “horizontal,  upi>er”,  “horizontal,  lower”,  “vertical, 
left”  and  “vertical,  right”,  respectively.  For  instance,  for  the  16  x  16  lattice,  the  sequences 
are  shown  in  Figure  15  and  defined  below:^^ 

0O,NW,hu(^)  =^(0,  fc), 

corresponding  to  the  points  denoted  with  Vi 

the  combination  of  V  and  <J , 

and  the  combination  of  v  and  t>  in  Fig.  15. 

(92) 

0O,NW,vr(k)  =  z(k,7), 

corresponding  to  the  points  denoted  with  <I , 

the  combination  of  <3  and  Vi 

and  the  combination  of  <3  and  A  in  Fig.  15. 

(93) 

lh,NW,hli^)  =«(7,A:), 

corresponding  to  the  points  denoted  with  A, 

the  combination  of  A  and  <3 , 

and  the  combination  of  A  and  t>  in  Fig.  15. 

(91) 

l^0,NW,vl{k)  =  z{k,0). 

corresponding  to  the  points  denoted  with  t>, 

the  combination  of  I>  and  V- 

and  the  combination  of  t>  and  A  in  Fig.  15. 

(95) 

will  use  s(i,j}  to  denote  the  value  of  the  MRF  at  lattice  site  (i,  j).  If  the  lattice  has  Tx  rows  and  T2 
columns,  then  (i,j)  €  {0, 1,  •  •  • ,  Ti  —  1}  x  {0, 1,  •  •  • ,  r2  —  1}. 
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Po,NE,hu(k) 

= 

z{0,k  +  8) 

(96) 

Po,NE,vrik) 

= 

z{k,  15) 

(97) 

Pa,NE,hl(k) 

= 

-(7,  A: +  8) 

(98) 

Po,NE,vl{k) 

= 

zik,8) 

(99) 

/3o,SE,huik) 

= 

^:(8,  A:  -|-  8 ) 

(100) 

l3o,SE,vr{k) 

=: 

z{k  4-  8, 15) 

(101) 

Po,SE,h0) 

= 

z(15,k  +  8) 

(102) 

Po,SE,vl{k) 

= 

z{k  -f-  8,8) 

(103) 

P(i,SW,hu{k) 

2^(8,  A:) 

(104) 

Po,SW,vrik) 

= 

~(A;  -f  8, 7) 

(105) 

l3o,sw,hi{k) 

= 

2(15,  A:) 

(106) 

0o,sw,vii  A" ) 

= 

z(^k  -(-  8, 0) 

(107) 

for  A;  =  0, 1,  •  •  • ,  7.  Note  there  is  overlap  in  the  sequences  For  instance,  l3o,NW,hu  and 

iio.NW,vt  both  contain  the  value  of  the  process  at  (0,0),  and  this  fact  is  reflected  in  Figure 
15  by  the  presence  of  both  V  and  t>  at  this  lattice  point. 

Let  us  now  consider  the  simplest  of  our  approximate  models.  Specifically,  we  take  as 
the  state  of  the  approximate  representation  just  the  averages  of  the  sequences  The 

state  at  any  node  then  has  sixteen  components: 


where: 


a;(s)  = 


•CivrvC-s) 
xne{s) 
xse{  -5) 


.Ti(.S) 


WoPs,i,hu 

lTo/?s,i,vr 

WoPs,i,hl 

lFo/^s,tVui 


(108) 


(109) 


for  i  £  {  NW,  NE,  SE,  SW}  and  where  WoPs,i,j  denotes  the  average  of  the  sequence  pa,i,j(k). 
Given  the  definition  of  the  state  (108), (109)  (which  will  be  generalized  shortly  to  allow 
general  wavelet  transform  approximations  to  the  sequence  (3s,i,j)',  the  conditional  parent- 
offspring  pdf’s  need  to  be  obtained  from  the  MRF  being  approximated.  Instead  of  using 
these  directly,  we  make  an  additional  approximation.  Lei  us  define  the  dowusliift  oper¬ 
ators  ai,  i  £  {NW,NE,SE,SW},  which  are  the  counterparts  of  the  upshift  operator  7 
defined  previously  (see  Figure  12).  In  particular,  we  denote  the  four  offspring  of  node  .s  as 
sai,i  £  {NW,NE,SE,SW},  where  the  subscript  refers  to  the  relative  spatial  location  of 
the  offspring  nodes.  In  the  exact,  non-redundant  representations,  the  following  relationship 
holds: 


Pz(t),ter,ai\4^hr£rsi^t^^  €  rscx,\Zr,r  er^)- 

Pz(t),terscci\z{T),Ter,,i{Zut  £  TaailZr^r  £  Fsy) 


(110) 
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for  i  e  {NW,NE,SE,SW}.  What  (110)  says  is  that  the  conditional  pdf  for  the  state  at 
node  .soi  depends  only  on  a  subset  of  the  values  making  up  the  state  at  the  parent  node  .s. 
For  example,  in  the  case  of  the  iVlF  offspring  of  node  s,  the  state  in  the  exact  representation 
at  node  samv  (that  is,  z(t),t  G  Fs^jw,,-)  depends  only  on  the  iVVF  component  of  the  state 
at  node  s  (that  is,  on  the  values  z(t),t  €  rs,ArM")-  Thus,  in  the  exact  representation  the 
state  at  node  samv  is  independent  of  the  values  of  the  MRF  at  the  points  in  Vs,ne^^s,se 
and  r5,.svv5  given  the  values  at  Irr  contrast,  it  is  not  true  in  general  in  the  sim¬ 

ple  approximate  representations  just  described  that  the  state  xlsamv)  is  independent  of 
xne(.?),  xse{^)  and  .rsvvls),  given  .'Cjvrv(s).  That  is,  simply  knowing  the  average  value  of 
a  process  along  each  side  of  a  square  region  does  not  completely  decorrelate  the  values  of 
the  field  inside  and  outside  the  region.  Nevertheless,  in  our  approximate  modehng  frame¬ 
work  we  will  make  exactly  this  assumption.  More  generally  and  precisely,  our  approximate 
modeling  methodology  yields  a  sequence  of  models  corresponding  to  differing  resolution 
approximations  to  the  boundary  processes  where  (108)  -  (109)  corresponds  to  the 

coarsest  of  these.  Using  the  same  symbols  Xi(s)  and  .r(s)  to  denote  the  state  components 
and  state  of  any  of  these  models,  we  construct  our  model  by  making  the  approximation 
corresponding  to  assuming  that  the  conditional  independence  property  holds,  i.e.  that: 

jPjl?(»ai)ti-(s)(^'ao;i  l-^s)  ?\'(sOi)|*j(s)(-Tsa;  l-^d-^))  (111) 

Since  the  field  being  approximated  is  assumed  to  be  jointly  Gaussian,  the  conditional 
density  function  (111)  is  parameterized  by  conditional  means  and  covariances  as  in  (.33)  - 
(35): 


Pa;(scv,;)|a;i («)(  -Tso;j  j..T}( ’5)  )  —  i  ^’scvi  ?  Fsa-,  )  (11^) 

where: 

Xsa.  =  E{,'i-(5ai)j.rj(s)} 

=  E{.'i:(saj).Ti(s)^}(E{,ri(s);i-i(.s)^})"'^A'd.s)  (113) 

Psai  —  E{(2:(sQ'j)  )(x’(sQ;,)  •'r’sQ';)  } 

=  E{a:(5Q'i);(;(5ai)^} 

-E{a;(sai).ri(s)^}(E{.ri(s).rj(s)^})~’^(E{.T(sa,)ri(s)^})^ 

(114) 

One  can  then  derive  the  matrices  A(s),i?(s)  and  Pq  in  the  multiscale  representation  of 
the  random  field: 


A(  .so'jvm )  = 

[A'l.  0,0.0] 

(115) 

A{  sone  )  = 

[0,A'2,0,0] 

(116) 

A{sasE)  = 

[0,0,A3,0] 

(117) 

/ll-scrsp/)  = 

[0,0,0,A4] 

(118) 

E{.T(sai)(a;i(.s))^ 

}(E{;f4s)(.Tj(.s))^})“^ 

(119) 

41 


Also: 


B(sai)B{8ai)^  —  Fsat  (120) 

Po  =  Efa’o-'i'o}  (121) 

The  assumption  (111)  is  directly  reflected  in  (115)  -  (118).  In  particular,  the  state  x{stti) 
is  a  function  only  of  the  component  of  the  parent  (cf.  (108)).  Thus,  the  assumption  in 
(111)  leads  to  relatively  simple  level-to-level  interpolations.  Indeed,  if  the  MRF  is  station¬ 
ary,  from  symmetry  we  see  that  not  only  do  the  parameters  /l(.s),i?(.s)  depend  only  on  the 
scale  of  node  s,  but  also,  A’l  =  K2  =  K's  —  Ka-  Thus,  in  this  case,  the  representations  are 
cpiite  simply  described,  and  more  importantly,  this  simple  structure,  in  addition  to  the  sub¬ 
stantially  reduced  dimensionality  of  the  approximate  representations,  leads  to  considerable 
efficiencies  for  smoothing  [17,  18,  19]  and  likelihood  calculation  algorithms  [43]. 

As  we  have  indicated,  the  generalization  of  the  coarsest  approximate  model,  with  state 
given  by  (108),  (109)  corresponds  to  using  wavelet  transforms  to  obtain  different  resolution 
representations  of  the  secpiences  l3s,i,j(k).  We  iitilize  the  wavelet  transform  for  discrete 
sequences  as  described  in  [6].  The  wavelet  transform  of  /3s,ij(k),  Ar  €  {1, 2,  •  •  • ,  is  a 

set  consisting  of  a  single  “scaling”  coefficient  and  —  1  “detail”  coefficients^^.  These 

are  computed  recursively  according  to*^: 

n=2Af 

fk  ~  X]  ^^^fn+2k-2  (122) 

n=l 

n=2M 

~  XI  37ifi+2k-2  (123) 

n— 1 

where  the  scaling  coefficients  and  detail  coefficients  are  and  respectively,  hn,gn  are 
impulse  responses  of  quadrature  mirror  filters  [26,  55]  of  length  2M,  and  where  _ 

/3s.i.j{k).  We  say  that  a  //^-order  representation  of  the  sequence  /3s^ij{k)  is  a  set  consisting 
of  the  scaling  coefficient  and  the  wavelet  coefficients  up  to  order  p  in  the  wavelet  expansion, 
and  that  a  zeroth-order  representation  is  a  set  consisting  of  just  the  scaling  coefficient.  We 
denote  the  operator  which  maps  the  sequence  l3s,i.j{k)  to  its  p‘*-order  representation  as 
Wp.  Note  that  \i p  =  N  —  m{s)  the  representation  is  complete,  since  it  contains  the  scaling 
coefficient  and  all  of  the  wavelet  coefficieirts.  For  p  >  N  —  7n{s)  we  take  Wp  =  l'FAr-m(s) 
(i.e.  if  there  are  fewer  than  p  scales  of  wavelet  coefficients,  we  keep  all  of  them). 

The  generalization  of  the  approximate  representation  based  on  averages  of  the  1-D 
sequences  I3g^i^j{k)  discussed  previously  now  just  involves  a  new  definition  for  the  state 

be  concrete,  we  assume  that  the  wavelet  transform  filter/do^vnsample  operations  are  iterated  until 
the  sequence  of  scaling  coefficients,  i.e.  the  downsampled  output  of  the  lowpass  component  of  the  wavelet 
filter  bank,  is  of  length  one.  More  generally,  one  could  stop  at  any  point  in  t  he  decomposition. 

*®Our  notation  is  slightly  different  from  that  in  [6].  In  particular,  in  [6],  increasing  superscript  j  cor¬ 
responds  to  lower  levels  in  the  decomposition  (i.e.,  fewer  wavelet  and  scaling  coefficients),  while  here  it 
corresponds  to  higher  levels. 
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variables  x{s).  In  particular,  simply  replace  (109)  with: 


Xi{s)  = 


^'^pl^s,i,vr 

^^pl^a,i,vl 


(124) 


where  14p/3s,i,j  denotes  the  p*^-order  representation  of  the  sequence  lis,i,j(k)  (a  vector  of 
length  2'^  if  /?  <  JV  -  771(5)  and  of  length  p  >  N  —  m{s)).  Thus,  the  state  at 

any  given  node  consists  of  sixteen  components,  each  a  27*'‘-order  representation  of  one  of 
the  1-D  boundary  sequences  Ps,i,j(k)  associated  with  the  state  a;(5).  Using  this  generalized 
definition  for  the  state,  and  making  the  assumption  in  (111),  the  parameters  ^(s),  B{s)  and 
Po  are  again  given  by  (115)  -  (121). 

Several  comments  are  in  order.  First,  note  that  a  simple  generalization  of  the  above  rep¬ 
resentation  would  be  to  allow  different  levels  of  approximation  for  different  components  of 
the  boundary  sequences  (e.g.  one  might  use  a  pj^-order  approximation  for  “vertical”  bound¬ 
ary  sequences  l3s,i,j,j  €  {ur,u/}  and  a  p^^-order  approximation  for  “horizontar  boundary 
sequences  lds,i,j^3  €  {/?.«, /i/}).  Examples  of  such  a  generalization  will  be  given  in  the  next 
section  in  the  context  of  approximate  representations  for  MRF  texture  models. 

Second,  note  that  even  if  all  of  the  wavelet  coefficients  are  retained  at  all  levels  (i.e. 
if  the  boundary  representations  are  complete),  the  representation  we  have  just  described 
will  be  exact  only  if  the  GMRF  is  Markov  with  respect  to  either  the  first  or  second-order 
neighborhood  in  Figure  11.  As  we  have  discussed,  higher-order  neighborhoods  lead  to 
thicker  boundaries,  and  this  leads  naturally  to  the  idea  of  taking  wavelet  expansions  of 
boundaries  of  width  two  or  more,  and  utilizing  these  as  the  state.  When  the  boundaries  are 
expanded  to  have  a  width  of  q  lattice  sites,  the  state  at  node  s  will  be  broken  up  into  the  p*^- 
order  representations  of  16?  sequences  of  length  With  this  expanded  family,  the 

approximate  representations  can  be  made  exact  for  any  GMRF  by  keeping  complete  wavelet 
expansioirs  of  all  boundary  sequences  (3s,i,i(k)  at  all  scales.  An  example  of  an  approximate 
representation  which  keeps  wavelet  coefficient  along  boundaries  of  width  two  is  discussed  in 
the  next  section. 

Third,  note  that  the  covariance  matrices  required  in  (119)  and  (120)  are  not  invertible  if 
the  representation  of  the  1-D  boundary  sequences  is  complete,  due  to  the  fact,  as  mentioned 
previously,  that  these  sequences  overlap.  For  instance,  in  Figure  15,  both  0Q,NW,hu 
0o,NW,vi  contain  the  value  of  the  state  at  pixel  location  (0, 0).  In  this  case  we  have  redundant 
information  and  hence  the  conditional  expectation  and  error  covariance  formulas  must  be 
modified  to  deal  with  this.  This  modification  is  a  straightforward  matter  as  discussed  in 
[50,  58]. 

Fourth,  note  that  not  only  has  the  dimensionality  of  the  representations  been  reduced 
in  going  from  the  exact  to  the  approximate  representations,  but  it  has,  in  fact,  been  made 
constant  at  the  first  N  —  p  levels  of  the  quadtree,  where  p  is  the  order  of  the  approximation 
and  the  MRF  is  defined  on  a  2^  X  2^^"  lattice.  In  particular,  the  dimension  of  the  state  at 
node  .?  is  equal  to  16  *  2^,  for  m{s)  <  N  —  p.  When  m  =  N  —  p,  the  boundary  sequence 
representations  are  complete  and  the  dimension  of  the  state  drops  by  a  factor  of  2  at  each 
level  thereafter. 
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Fifth,  because  our  approximate  models  keep  only  limited  resolution  versions  of  MRF 
values  along  1-D  boundaries,  the  quadtrees  for  these  approximate  models  may  require  more 
levels  than  the  exact  model.  For  example,  consider  an  MRF  over  a  4  x  4  region.  The 
exact  representation  of  this  field  in  our  framework  has  only  a  single  level,  since  the  exterior 
boundaries  and  mid-lines  form  the  entire  4x4  region  (consider  Figure  15  adapted  to  a  4  x  4 
grid).  On  the  other  hand,  a  first-order  Haar  approximation  would  retain  only  the  sixteen 
average  values  of  pairs  of  horizontal  or  vertical  points  at  the  first  level,  only  twelve  of  which 
are  independent  thanks  to  the  overlap  in  the  sequences.  Consequently,  in  this  case,  we 
need  a  second  level,  corresponding  to  “averages”  of  single  points,  to  completely  represent 
the  field. 

Finally,  the  order  of  the  approximations  required  to  achieve  a  desired  level  of  fidelity 
in  the  approximate  model  depends,  of  course,  on  the  statistical  structure  of  the  specific 
GMRF  under  study.  In  the  next  section  we  present  examples  which  illustrate  this  for 
several  GMRF’s  and  a  number  of  different  approximate  representations. 

5  Examples  of  Approximate  2-D  Gaussian  MRF  Represen¬ 
tations 

Ill  this  section  we  illustrate  sample  paths  generated  by  our  approximate  representations  for 
two  examples  of  separable  Gaussian  MRF’s  and  then  for  two  examples  of  non-separable 
Gaussian  MRF’s.  Separable  MRF’s  were  one  of  the  first  widely  used  image  models  because 
of  their  simple  covariance  structure  [34],  while  non-separable  GMRF’s  have  been  widely 
used  in  the  context  of  texture  representation  [13,  14,  22,  23,  24,  47,  48]. 

5.1  Separable  Gaussian  MRF  Examples 

Consider  a  separable  Gaussian  MRF  defined  on  a  2^  x  2^  lattice  with  a  covariance  function 
given  by: 

E{zii,j}zik,l)}  =  (i25) 

where  and  py  are  one-step  correlation  parameters  in  the  vertical  and  horizontal  directions. 
These  random  fields  are  Markovian  with  respect  to  the  second-order  neighborhood  given  in 
Figure  11. 

Let  us  construct  a  zeroth-order  Haar  approximate  representation  of  this  MRF.  In  this 
case,  the  state  at  the  root  node  of  the  tree  consists  of  sixteen  values,  representing  the 
averages  across  each  of  the  1-D  components  of  the  boundary.  Since  this  state  variable  is 
just  a  linear  function  of  the  values  of  the  random  field,  its  covariance  sttncture  can  l)e 
calculated  directly  from  knowledge  of  the  averages  taken  and  the  covariance  function  in 
(125).  Indeed,  the  covariance  matrix  for  the  four  averages  corresponding  to  the  north-west 
corner  of  the  state  at  the  root  node  is  given  by; 

Wof3o,NW,hu  Wo/3o^NW,hu  0 

IFo/lOiAflVj'ur  II' ^  V  Qy'^x 

WQl5o^mv,hi  Wol3Q^mv,hi  Qx>^y  i'  ^y  t 
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where: 


a;,  =  {2^-‘(l+/7,)/{l-^,)-2^Vl-/,f“')/(l-/,,)2  (127) 

u;,:  =  (2^''-i(  1  +  )/( 1  -p,:)-  2p.A  1  -  pf~' )/( 1  -  p^f  ( 128 ) 

=  {(I-  pf -')!({-  p,]]((l-  pf~^)l(l-  p^))  (129) 

Qr.  =  pI  '  (130) 

Qy  =  Py  (131) 

Figures  17a  ajul  17b  illustrate  256  X  256  sample  paths  of  exact  representations  of  the  sep¬ 
arable  randoni  fields  for  p^;  —  py  =  0.7  and  p^.  =  py  —  0.9,  respectively.  Sample  paths 
of  zeroth-order  Haar  approximations  of  these  MRF's  are  shown  in  Figures  17c  and  17d, 
respectively Note  that  for  p^  —  Py  =  0.7,  the  zeroth-order  approximation  is  visually 
similar  to  the  exact  representation,  with  only  minor  boundary  effects  apparent,  caused  by 
the  approximation  in  the  first-order  Haar  model,  i.e.  the  neglecting  of  the  residual  correla¬ 
tion  in  adjoinijig  regions  when  we  are  given  only  the  coarse  Haar  approximation  of  the  field 
along  common  boundaries.  As  the  coupling  between  pixels  is  increased,  the  effects  of  this 
coarse  approximation  become  more  apparent,  as  seen  in  the  first-order  approximation  of  the 
.separable  MHF  with  p^  =  py  =  0.9.  This  indicates  that  such  fields  wifi  in  general  require 
higher-order  approximations.  We  defer  the  illustration  of  such  higher-order  approximations 
of  fields  to  the  following  subsection,  in  which  we  describe  several  examples  of  the  use  of  our 
modeling  methodology  to  represent  natural  textures. 


5.2  Non-Separarable  Gaussian  MRF  Examples 


(  'onsider  the  class  of  GMRF’s  defined  by  the  following  2-D  autoregressive  model  [12,  36]: 

H  nMi-  kj  +  (132) 

(W)€D 


where  ri.,i  =  is  <*  neighborhood  around  the  origin  (0,0),  the  Gaussian  driving  noise 

p(uj)  is  a  locally  correlated  sequence  of  random  variables,  and  (i,j]  €  {0, 1,  •  •  • , 2^^^}  X 
{0. 1,  •  •  • ,  2^}.  In  the  examples  below,  the  set  D  corresponds  to  the  neighborhood  .sets 
in  Figure  11.  For  instance,  the  set  corresponding  to  a  first-order  neighborhood  is  Z)  = 
{(0, 1),  (0,  — 1),(1,0),(— -1, 0)}.  In  addition,  we  interpret  the  2^  X  2^  lattice  as  a  toroid,  i.e. 
the  independent  variables  {i,j)  in  (132)  are  interpreted  modulo  2^.  For  instance,  the  first- 
order  neighborhood  of  lattice  site  (0,0)  is  given  by  the  set  {(1,0),  (0, 1),  (0,2^  —  1),  (2^  — 
1,0)}.  Finally,  the  correlation  structure  of  the  driving  noise  is  given  by: 


E{e{i,j)€(i-  k,j  -  /)} 

and  has  the  property  that: 


E{e(f,i)r(A:,/)}  = 


(  if  A-  =  /  =  0 

<  a,/  if(A-./)eii 

[  0  if(A'./)^  .D 

(T^  for  /  =  k.j  =  / 

0  else 


(133) 


(134) 


'■'We  emphasize  that  the  Figures  ]7c  -  17d  depict  sample  paths  oi  the  approximate  representations,  and 
not  approximate  representations  of  the  sample  paths  in  Figure  17a  -  17b. 
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Figure  17:  Sample  paths  of  a  separable  MRF  with  correlation  structure  'E{z{  i,j)z{k,l)}  = 
Pj:  —  Py  —  0-7  and  px  =  Py  ~  0.9  are  shown  in  (a)  and  (b)  respectively. 
Sample  paths  of  zeroth-order  approximate  representations  of  these  fields,  based  on  the  Haar 
wavelet,  are  shown  in  (c)  and  (d).  The  stronger  correlation  between  neighboring  pixels  in 
(b)  leads  to  boundary  effects  in  the  sample  path  of  the  approximate  representation  shown 
in  (d),  which  can  be  eliminated  by  using  higher-order  approximations. 


From  (134),  and  the  fact  that  the  random  field  is  Gaussian,  one  can  prove  that  the  autore¬ 
gressive  model  above  does  imply  that  z{i,j  )  is  a  Mai'kov  Random  Field  [64],  We  refer  to  the 
model  (132)  as  a  Q*^-order  MRF  if  the  set  D  corresponds  to  the  Q^^'-order  neighborhood 
of  Figure  11. 

Infinite  lattice  versions  of  the  processes  in  ( 132)  were  introduced  in  [64]  and  their  toroidal 
lattice  counterparts  have  been  thoroughly  studied  in  the  context  of  texture  representation 
[13,  14,  22,  23,  24,  47,  48].  The  correlation  structure  of  these  MRF’s  cannot  be  explicitly 
written  down  as  in  the  previous  example.  However,  the  specific  statistics  and  correlations  (as 
in  ( ll.l)  -  (118))  requii-ed  to  construct  our  multiscale  approximate  models  can  be  computed 
efficiently  using  2-D  FFT’s  because  of  the  fact  that  correlation  matrices  for  these  random 
fields,  assuming  lexicographic  orderiirg,  are  block  circulant  with  circulant  blocks  and  hence 
these  random  fields  are  whitened  by  the  2-D  Fourier  transform  [36].  In  particular,  denote 
by  z  the  set  of  values  z(t)  stacked  into  a  vector  (with  lexicographic  ordering),  and  denote 
the  correlation  matrix  of  the  MRF  by  R^z-  Then: 

FRzzF*  =  A  (135) 

where  F  is  the  2-D  Fourier  transform  matrix  and  A  is  a.  diagonal  matrix  of  the  eigenvalties 
of  Rzz-  The  wavelet  coefficients  required  in  the  approximate  representation  correspond  to 
linear  functions  of  the  values  z(t).  That  is: 

^Fpf3s,i,j  =  WpSz  (136) 

=  Lz  (137) 

where  the  matrix  L  is  the  product  of  the  wavelet  transform  operator  Wp  and  a  “selection” 
matrix  S  which  generates  the  vector  from  z.  Thus,  to  compute  the  correlation  matrices 
required  in  the  approximate  representation,  we  need  only  compute  functions  of  the  form: 

LiRzzL^  =  (LiF*)A(FL2)  (138) 

Indeed,  as  described  in  [4.3],  the  structure  of  the  approximate  representations  and  the  sta- 

tiouarity  of  the  GMRF  allow  us  to  compute  the  required  correlations  with  only  2^  2-D 
Fourier  transform  operations  per  level  of  the  representation,  where  p  is  the  order  of  the 
approximation.  Furthermore,  these  calculations  need  only  be  performed  once,  since  they 
are  used  simply  to  determine  the  parameters  in  the  multiscale  approximate  model. 

Figure  18a  illustrates  the  “wool”  texture  from  [14].  Three  sample  paths  of  approximate 
representations  of  this  field  based  on  the  Haar  wavelet  are  shown  in  Figures  18b  to  18d. 
The  wool  texture  is  cori’esponds  to  a  fourth-order  version  of  the  model  (132),  with  the 
coefficients  given  in  Table  1. 

Figures  18b  and  18c  correspond  to  zeroth  and  first-order  approximations,  respectively. 
Note  in  Figure  18c  that  some  of  the  boundary  effects  apparent  in  Figure  18b  have  disap¬ 
peared,  due  to  the  increase  in  the  approximation  order.  Since  the  wool  texture  is  actually 
Markov  with  respect  to  the  fourth-order  neighborhood  of  Figure  11,  an  exact  representa¬ 
tion  of  this  field  woidd  reciuire  that  boundaries  of  width  equal  to  two  be  kept.  To  this  end, 
the  approximation  shown  in  Figure  18d  takes  as  state  variables  first-order  approximations 
to  these  double  width  boundaries.  Essentially  all  of  the  boundary  effects  present  in  the 
Figures  18b  to  18c  have  been  eliminated,  and  this  representation  appears  to  have  retained 
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Figure  18:  A  sample  path  of  a  Gaussian  MRF  representing  the  “wool”  texture  of  [14]  is 
shown  in  (a).  Figures  18b  -  ISd  illustrate  sample  paths  of  approximate  representations  of 
this  MRF  based  on  the  Haar  wavelet.  Zeroth  and  first-order  approximations  are  used  in 
(b)  and  (c), respectively.  In  (d),  a  first-order  approximation  based  on  boundaries  of  width 
two  is  used.  Note  that  as  the  order  of  the  approximate  model  is  increased,  the  boundary 
effects  disappear,  and  that  for  relatively  low-order  models  an  approximate  representation 
which  retains  most  of  the  qualitative  and  statistical  features  of  the  original  MRF  can  be 
obtained. 
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{k,l) 

hk,l 

(k,l) 

hk,i 

(1,0) 

0.4341 

(0,2) 

0.0.592 

(0,1) 

0.2182 

(-1,2) 

-0.0302 

(-1,1) 

-0.0980 

(1,2) 

-0.0407 

(1,1) 

-0.0006 

(-2,1) 

0.0406 

(2,0) 

-0.08.36 

(2,1) 

-0.0001 

Table  1:  Coefficients  in  the  model  (132)  for  the  “wool”  texture. 


ikj) 

hk,t 

(k,l) 

hk,l 

(i,0) 

0.5508 

(0,2) 

0.0139 

(0,1) 

0.2498 

(-1,2) 

-0.0085 

(-1,1) 

-0.1164 

(1,2) 

-0.0058 

(1,1) 

-0.1405 

(-2,1) 

-0.0008 

(2,0) 

-0.0517 

(2,1) 

0.0091 

Table  2:  Coefficients  in  the  model  (132)  for  the  “wood”  texture. 


the  essential  statistical  and  qualitative  features  of  the  exact  representation  used  to  generate 
Figure  18a. 

Figure  19a  illustrates  the  “wood”  texture  from  [14],  and  three  approximations  of  this 
MRF  based  on  the  Haar  wavelet  are  shown  in  Figures  19b  -  19d.  The  wood  texture  corre¬ 
sponds  to  a  fourth-order  version  of  the  model  (132),  with  the  coefficients  given  in  Table  2. 

This  texture  clearly  has  a  very  asymmetric  correlation  structure,  and  thus  we  represent 
the  vertical  and  horizontal  boundary  with  different  levels  of  approximation.  In  Figure  19b, 
the  horizontal  and  vertical  boundaries  are  represented  with  second  and  zeroth-order  approx¬ 
imations  respectively.  In  Figures  19c  and  19d,  the  horizontal  boundaries  are  represented 
with  fourth  and  sixth-order  approximations,  respectively,  whereas  the  vertical  Iroundary  is 
still  represented  with  a  zeroth-order  approximation.  As  the  complexity  of  the  representation 
increases,  the  sample  paths  of  the  approximate  random  fields  have  fewer  boundary  effects. 
The  approximate  representations  used  to  generate  Figures  19c  and  19d  appear  to  accurately 
represent  the  qualitative  and  statistical  features  of  the  MRF.  An  interesting  point  here  is 
that  the  level  of  representation  only  needs  to  be  Increased  in  one  direction  to  obtain  an  ex¬ 
cellent  representation  of  the  field.  Also,  the  neighborhood  of  this  MRF  is  fourth-order  (see 
Figure  11)  and  thus  double  width  boundaries  would  be  needed  in  an  exact  representation. 
The  fields  shown  in  Figures  19b  to  19d,  however,  use  only  the  thinner  boundaries  in  forming 
states.  Several  experiments  were  done  in  which  we  used  the  double  width  Ijoundaries  in 
forming  states  for  models  analogous  to  those  in  Figures  19b  to  19d.  It  was  found,  how¬ 
ever,  that  there  were  no  visual  differences  between  the  single  and  double  width  approximate 
representations. 

Three  approximations  of  the  “wood”  texture  based  on  the  Daubechies  8  wavelet  de- 
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Figure  19:  A  sample  path  of  a  Gaussian  MRF  representing  the  “wood”  texture  of  [14]  is 
shown  in  (a).  Figures  19b  -  19d  illustrate  sample  paths  of  approximate  representations  of  the 
MBF  based  on  the  Haar  wavelet.  The  structure  of  the  MRF  suggests  using  approximations 
which  use  relatively  low  order  representations  of  vertical  Imundaries.  The  approximate 
representations  used  to  generate  Figures  19b  -  19d  used  zeroth-order  representations  of  the 
vertical  boundaries,  and  second,  fourth  and  sixth-order  representations  for  the  horizontal 
boundaries,  respectively. 
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Figure  20:  Figures  20a  -  20c  illustrate  sample  paths  of  approximate  representations  based 
on  the  Daubechies  8  wavelet,  with  the  same  orders  of  approximation  as  in  Figures  19b  - 
19d. 

scribed  in  [26]  are  illustrated  in  Figures  20a  to  20c.  The  order  of  the  approximations  are 
identical  to  the  orders  for  Haar  approximations  in  the  previous  example.  We  note  that 
there  is  no  apparent  difference  between  the  approximations  based  on  the  Haar  wavelet  and 
the  Daubechies  8  wavelet.  That  is,  at  least  for  this  example,  and  for  the  others  we  have 
examined,  the  critical  issue  in  model  fidelity  appears  to  be  model  order  rather  than  the  par¬ 
ticular  choice  of  the  wavelet  used.  Furthermore,  as  these  examples  indicate,  we  can  achieve 
quite  high  quality  results  with  low-order  models,  which  in  turn  lead  to  extremely  efficient 
algorithms  as  in  [16,  17,  18,  19].  In  addition,  as  we  briefly  discuss  in  the  next  section, 
for  GMRF’s  which  have  particular  directions  in  which  correlation  structures  are  oscillatory 
rather  than  monotonically  decaying  (such  as  the  one  describing  the  “wood”  texture),  there 
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may  be  different  choices  of  bases  other  than  wavelets  that  lead  to  high  fidelity  models  of 
even  lower  dimension. 

6  Discussion  and  Conclusions 

In  this  paper,  we  have  shown  how  to  represent  reciprocal  and  Markov  processes  in  one 
dimension  and  Markov  random  fields  in  two  dimensions  with  a  class  of  multiscale  stochastic 
models.  This  modeling  structure  provides  a  framework  for  the  development  of  efficient, 
scale-recursive  algorithms  for  a  variety  of  statistical  signal  processing  problems.  The  rep¬ 
resentations  in  1-D  rely  on  a  generalization  of  the  mid-point  deflection  construction  of 
Brownian  motion.  In  2-D,  we  introduced  a  “mid-line”  construction  which  leads  to  a  class 
of  models  with  scale-varying  dimension.  Since  for  reasonable  size  fields  this  state  dimen¬ 
sion  may  be  prohibitively  large,  we  also  introduced  a  class  of  multiscale  approximate  MllF 
representations  based  on  1-D  wavelet  transforms  of  the  MRF  along  1-D  boundaries  of  mul¬ 
tiresolution  partitionings  of  the  2-D  domain  of  interest.  This  family  of  models  allows  one  to 
tradeoff  complexity  and  accuracy  of  the  representations,  and  i>rovides  a  framework  for  the 
development  of  extremely  efficient  estimation  and  likelihood  calculation  algorithms.  Ex¬ 
amples  demonstrated  that  for  relatively  low-order  models,  an  approximate  representation 
which  retains  most  of  the  qualitative  and  statistical  features  of  the  original  MRF  can  be 
obtained.  Moreover,  these  approximate  models  lead  to  algorithms  which  have  a  complexity 
that  is  constant  per  pixel,  in  contrast  to  the  algorithms  normally  associated  with  MRF’s 
which  have  a  per  pixel  complexity  which  increases  with  image  size. 

We  feel  that  the  results  presented  In  the  preceding  section,  together  with  the  substantial 
flexibility  of  the  multiscale  modeling  framework,  both  demonstrate  the  promise  of  this 
framework  for  image  and  multidimensional  signal  processing  and  also  suggest  a  rich  set 
of  questions  for  further  investigation.  F'or  example,  our  comparisons  here  between  exact 
and  approximate  representations  have  focused  only  on  their  visual  characteristics.  It  is 
important  to  develop  other  measures  of  the  quahty  of  these  fields  and  in  particular  to 
define  and  analyze  metrics  that  allow  us  to  place  the  trade-off  between  fidelity  and  model 
order  on  a  rational  basis.  For  instance  one  natural  basis  for  measuring  the  distance  between 
two  models  corresponds  to  analyziirg  error  probabilities  associated  with  a  binary  hypothesis 
test  in  which  one  must  decide  whether  a  given  random  field  is  a  sample  path  of  the  exact 
or  approximate  representation.  In  particular,  if  the  probability  of  error  is  near  to  1/2,  then 
the  fields  are  nearly  indistinguishable. 

It  is  also  important  to  investigate  metrics  directly  related  to  performance  in  applica¬ 
tions  in  which  performance  based  on  MRF  models  and  on  approximate  models  of  varying 
orders  can  be  compared.  In  particular,  the  GMRF  texture  models  are  often  used  in  region 
labelling  and  image  segmentation  applications.  It  would  be  of  interest  to  develop  algo¬ 
rithms  for  these  applications  based  on  the  approximate  representations,  and  compare  them 
in  terms  of  performance  and  required  computation  to  the  standard  algorithms.  Note  in 
particular,  that  while  our  models  are  approximations  to  MRF's.  they  may  very  well  lead 
to  algorithms,  which  are  not  only  far  simpler  computationally  but  which  also  have  equal  or 
better  performance.  Specifically,  an  extremely  important  point  is  that  all  models,  including 
MRF’s  and  multiscale  models,  are  idealizations,  and  thus  it  is  not  at  all  obvious  that  MRF 
based  algorithms  will  perform  better  than  those  based  on  our  models.  Indeed,  given  the  fact 
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that  the  miiltiscale  modeling  framework  leads  to  efFicieiit  algorithms  and  is  richer  than  the 
class  of  MRF’s  (in  that  any  MRF  can  be  approximated  to  any  desired  degree  of  accuracy 
using  multiscale  models),  there  is  a  strong  argument  in  favor  of  the  framework  described 
here. 

In  addition,  there  are  numerous  other  directions  in  which  the  ideas  in  this  paper  can 
be  extended  and  further  developed.  For  example  there  are  strong  motivations  for  the  con¬ 
sideration  of  alternatives  to  wavelet  transforms  for  the  approximate  rei^resentations  used 
in  our  multiscale  models.  For  instance,  it  is  certainly  possible  to  consider  non-orthogonal 
multiresolution  approximations  to  the  values  of  MRF’s  along  1-D  boundaries.  Indeed,  one 
possibility  is  to  use  linear  or  higher-order  polynomial  interpolation  —  i.e.  to  perform  gen¬ 
eralized  mid-point  deflection  along  each  of  the  1-D  boundaries.  Also,  and  perhaps  more 
importantly,  wavelet  packet  basis  functions  [62]  may  provide  lower- dimensional  approxima¬ 
tions  for  some  GMRF’s,  such  as  the  “wood”  texture  discussed  in  the  previous  section,  in 
which  the  random  field  has  bandpass  (i.e.  oscillatory)  rather  than  low-pass  characteristics  in 
one  or  more  directions.  Indeed,  techniques  such  as  in  [25],  suggest  the  idea  of  choosing  the 
“optimal”  set  of  basis  functions  within  some  class,  where  optimality  might  be  measured  in 
terms  of  the  compactness  of  the  representation,  i.e.  in  terms  of  the  order  of  the  approximate 
model  required  to  achieve  an  acceptable  representation. 
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